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Automated Network Analyzers for the 
0.9- to 12.4-GHz Range 


By J. G. EVANS, F. W. KERFOOT, and R. L. NICHOLS 


(Manuscript received December 22, 1975) 


Two highly accurate automated network analyzers have been developed 
to cover the frequency ranges of 0.9 to 4.2 GHz and 8.7 to 12.4 GHz. In 
addition to the directly measured quantities of insertion loss and phase, 
the analyzers provide insertion delay, impedance, and two-port parameters 
as derived quantities. 

The functional architecture of the systems ts described along with some 
of the novel techniques used to obtain the high signal-to-noise ratios and 
system stabilities needed for the precise measurements of loss, phase, and 
delay. General descriptions of these techniques are given along with some 
of the design considerations. 

Experimental results that demonstrate measurement accuracy are pre- 
sented for one-port and two-port devices using both coaxial and waveguide 
connectors. 


l. INTRODUCTION 


New automated measuring instrumentation has been developed for 
making network and component characterizations in the microwave 
frequency range from 0.9 to 12.4 GHz. The development was stimu- 
lated by characterization requirements growing out of the design of 
new communication systems and the redesign of older systems. These 
systems are achieving greater capacity due to advances in technology, 
both in the areas of improved analysis and improved components. 
These advances have led to requirements on measurements at micro- 
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wave frequencies that heretofore had been needed only at lower 
frequencies. 

This new instrumentation was implemented as an extension to the 
Computer Operated Transmission Measuring Sets (cotms)! operating 
in the Transmission Technology Laboratories at Bell Laboratories. 
This implementation avoids the duplication of costly common equip- 
ment and allows the utilization of many of the proven COTMS measure- 
ment techniques. 

The microwave instrumentation, in concert with the coTms system, 
is intended for the characterization of linear two-port networks, either 
active or passive. Measurement terminals are provided for networks 
with either coaxial or waveguide terminals. Provision is made for 
complete four-parameter characterization of network elements such 
as transistors (e.g., H parameters, Y parameters). 

The full 0.9- to 12.4-GHz frequency range is covered by two micro- 
wave systems of very similar design. The lower frequency system, com- 
pleted in November 1972, covers the range 0.9 to 4.2 GHz.° Its first 
applications were in characterizing solid state components. The 
second instrument, completed in May 1971, covers the range 3.7 to 
12.4 GHz," spanning the most widely used common-carrier frequency 
bands. It found immediate use in characterizing filter elements, and 
later complete networks for the TN1 radio system. 

This paper reviews the measurement capabilities and the overall 
design of both measuring systems. Several system-design considera- 
tions of special interest and error mechanisms that limit accuracy are 
described. This paper concludes by presenting experimental results 
that illustrate the applications of the measurement systems and demon- 
strate their high measurement accuracy. 


ll. MEASUREMENT SET CHARACTERISTICS 
2.1 General 


The two microwave measuring sets which cover the frequency 
ranges of 0.9 to 4.2 GHz and 3.7 to 12.4 GHz are fully automated 
network analyzers that were specifically designed for high-accuracy 
measurements. The sets are extensions of the 50-Hz to 1-GHz Com- 
puter Operated Transmission Measuring Set (corms). The CoTMS sys- 
tem provides input frequencies for the microwave signal source, IF 
detection circuitry, and the control computer. Extensive use is made of 
the computation capability to enhance measurement convenience and 
accuracy. Measurements are made under computer direction and the 


“For brevity this system will be referred to as the 1- to 4-GHz measuring set. 
t For brevity this system will be referred to as the 4- to 12-GHz measuring set. 
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calculations for derived parameters are made in real time. Accuracy 
is enhanced by the correction of known systematic errors and the 
averaging of random errors. The resulting measurement speed and 
accuracy are significantly greater than those previously achieved with 
manual techniques. 


2.2 Measurements 


The basic quantities measured by the test sets are loss (or gain) 
and phase. Insertion loss and insertion phase of a network are deter- 
mined by subtracting the loss and phase obtained with the test set 
connectors directly connected from the loss and phase measured 
with the network inserted between the connectors. Special care was 
taken to achieve an uncertainty in test-port impedance of less than 
2 percent on both the coaxial and waveguide test ports in order to 
minimize impedance mismatch errors. 

Envelope delay, 7, is a derived parameter approximated from in- 
sertion-phase measurements at two closely spaced frequencies. 


— 06 be — 
* 3600F 360(F2 — Fi) 





for 6 in degrees. (1) 


The return loss or reflection coefficient of a single-port device under 
test (DuT) is also a derived parameter.? The measurement is made with 
the hardware switched (for example, see the S1; and S22 configurations 
in Fig. 4) to couplers on the coaxial test ports and magic tee hybrids 
on the waveguide test ports. Three measurements with independent 
reflection calibration standards are used to substantially reduce errors 
associated with impedance mismatches and poor coupler directivity.’ 

The four-parameter characterization of a puT entails a complicated 
set of calibration measurements and switching of the hardware into 
four independent measurement configurations. The four measure- 
ments are closely related to the scattering (S) parameters? and are 
easily transformed to other parameter representations. Determination 
of the S parameters entails the measurement of three reflection co- 
efficient standards on each test port and a transmission standard 
connecting the test ports. Automatic switching provides the proper 
interconnections for two transmission and two reflection measure- 
ments corresponding to S21, S12, S11, and S22. The measurements 
are followed by a calculation to convert S parameters that are refer- 
enced to test set impedances to S parameters referenced entirely to 
the calibration standards, thereby eliminating any impedance-mis- 
match errors associated with imperfect test set impedances.’ More 
details on the hardware used for the above measurements are pre- 
sented in later sections. 
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2.3 Test set performance 


The uncertainties of measurements made on both test sets are de- 
pendent upon random errors, uncorrected systematic test set errors, 
and uncertainties in the calibration standards. The random errors can 
be reduced by signal averaging. For a typical measurement, the user 
can specify an upper bound on the numerically computed variance of 
the mean of every directly measured quantity. Averaging slows down 
the measurement rate. The measurement rate is five frequency points 
per second taking a single observation. The rate decreases to one fre- 
quency point per second if 32 observations are averaged. The typical 
standard deviation of a single observation at low insertion loss is 
0.0017 dB and 0.01 degree. 

Table I gives a representative list of the quantities measured and 
their respective resolutions. The resolution, or differential accuracy, 
is important in the characterizations of the small loss and phase varia- 
tions in the passband of a filter. 

Systematic test set errors may give rise to uncertainties in absolute 
characterizations that are larger than the differential uncertainties. 
Some of these errors, such as detection nonlinearity, are functions of 
many parameters and are not practical to remove by computation. 


Table |— Test set performance 
Characteristic Range Resolution 
Frequency 0.9 to 4.2 GHz; 3.7 to 12.4 GHz in two | 1 Hz 
test sets. 
Test rate Five measurements per second with no 


averaging for insertion loss and phase 
Insertion measurements” 


Loss 0 to 80 dB 0.002 dB 
Phase 0° to 360° 0.01° 
Delay (0 to 1980)/AF nst 0.06/AF ns 
Reflection measurements” 
Reflection coefficient 0 tol 0.02% 
Two-port measurements” 
Reflection 
parameters (S11, S22) 0 tol 0.02% 
Transmission 
parameters (S12, S21) 0 to 80 dB 0.002 dB 
0° to 360° 0.02° 
Incident power level —15 to —55 dBm (0.9 to 4.2 GHz) 1.0 dB 


—5 to —30 dBm (3.7 to 12.4 GHz) 2.0 dB 


*Loss and phase random uncertainties increase as loss increases, decreasing the 
signal-to-noise ratio. The variance of a single measurement increases as follows: 


Floss = 0.0017 X 104/29 dB 
Tphase = 0.01 X 102/20 deg. 


Lis the loss in dB of the put. The random uncertainty can be reduced by averaging. 
t AF is the frequency difference in eq. (1) expressed in MHz. 
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They have been reduced by careful system design to a magnitude 
typically less than the relationships given below when the best avail- 
able calibration standards are used.* 


[|ASj7| < 0.0015(1 + F/12)|1 — Sjz? — Sij-Sjt| (2) 
ASij 
S17 





IIA 





0.0001(1 + F) + 0.0001 | Siz|— 
+ KL + F/12)({S11] + [S22[) (3) 


kK = 0.01 for a simple insertion measurement 
K = 0.002 for mismatch errors after computer correction. 


The parameter (/) is the measurement frequency expressed in giga- 
hertz. The family of parameters (S77) are the scattering parameters 
of the device under test. 

Equation (2) can be applied to a one-port reflection coefficient 
characterization by setting |S12-S21| equal to zero. Uncertainty in 
the characterization of the standards is the dominant contributor to 
this partially empirical relationship. 

Equation (3) is composed of a crosstalk uncertainty, an impedance 
mismatch uncertainty, and a residual uncertainty from many other 
mechanisms. When a simple insertion measurement is made without 
compensation for the deviation of the test set impedances from 
nominal values, then the coefficient on the last term must be increased 
from 0.002 to 0.01. A more comprehensive discussion of the error 
mechanisms is given in Section V. 


Ill. BLOCK DIAGRAM DESCRIPTION 
3.1 Overall description 


A simplified block diagram illustrating the relationship between 
cots and the broadband microwave test sets is shown in Fig. 1. The 
test sets are quite similar in system design. They consist of source, 
comparison unit, heterodyne converter, and interface circuitry. 

The cotTms source supplies a computer-controlled low-frequency 
input signal to source circuitry in each test set. The microwave mea- 
surement frequency is developed in the source circuitry by heterodyn- 
ing and frequency multiplication of the input signal. Signal purity 
of the output signal is enhanced by computer-controlled yie filters. 

The comparison circuitry interconnects the source and detector via 
a switching network that rapidly interchanges an unknown path, 
containing the put, and a standard path. Full two-port characteri- 


* Equation (2) is valid if the uncertainty in the reflection coefficients of the im- 
pedance standards is less than 0.001. Equation (3) is valid if the uncertainty in 
S12 and S21 of the transmission standard is less than 0.0001. 
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zation of the unknown is provided by additional switching circuitry 
in the unknown path. 

The heterodyne circuit provides a fixed intermediate frequency 
input signal to the computer-controlled corms detector. This signal 
contains information on the differences in amplitude and phase be- 
tween the unknown and standard paths. The heterodyne circuit also 
provides a phase-reference input signal to the corms detector at the 
same fixed intermediate frequency. 

The interface provides digital control of the test system circuitry 
by the corms computer and informs the computer of circuit status. 
All digital lines between cots and the test systems pass through the 
interface. Manual control of each line is also provided by a toggle 
switch for maintenance and other special purposes. 


3.2 Source 


The test frequency in each applique is developed by heterodyning 
a signal at frequency f; received from coTms with a signal at frequency 
fe (see Fig. 2). The signal f2 is derived from reference signals coherent 
with the synthesizer-derived signal f:. The signal at frequency f: + fe 
is applied to the input of a comb generator which produces a harmonic 


DEVICE 
UNDER 
TEST 





COMPARISON UNIT 










COTMS 
SOURCE 





MICROWAVE 
SOURCE 


HETERODYNE 


CONVERTER 
INTERFACE 
COTMS 
COMPUTER 


Fig. 1—Simplified block diagram of the coTms-microwave measuring systems. 
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4—12 GHz MEASURING SET : fp = 825 MHz, 75 MHz <f, < 300 MHz 


Fig. 2—Block diagram of the source circuitry. 


output spectrum extending beyond 13 GHz. A single sinusoidal har- 
monic is selected by four stages of computer-controlled yie filtering, 
producing a test signal at frequency F = N(f1+ fe). This weak 
signal is amplified with solid state amplifiers in the 1- to 4-GHz test 
set and with a traveling-wave-tube amplifier in the 4- to 12-GHz test 
set. 


3.3 Heterodyne section 


In the heterodyne section (see Fig. 3), the microwave source signal 
(F) is split to excite the put, generate a reference-phase signal, and 
generate the local-oscillator signal. The local-oscillator (Lo) signal 
and the signals at frequency F are combined to form fixed intermediate- 
frequency signals (1F) which are sent to the corms detection circuitry. 
These 1F signals are processed by coTMs to extract phase and ampli- 
tude information. 

The source signal applied to the upper path shown in Fig. 3 excites 
the put in the comparison circuit. The put modulates the amplitude 
and phase of this signal, which are ultimately detected in corms. To 
be compatible with the corms frequency range, it is necessary to 
heterodyne this signal with an Lo signal to generate a lower-frequency 
IF signal. 

A reference signal is necessary to make a phase measurement. This 
reference signal is generated in the lower path in Fig. 3 in an almost 
identical fashion as that described above. 

The Lo signal common to both 1F mixers is offset from the measure- 
ment frequency by 0.495 GHz. This signal is generated by the mixer 
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Fig. 3—Block diagram of the heterodyne circuitry. 


circuit shown in the middle path of Fig. 3. The fixed 0.495-GHz signal 
is derived from reference signals coherent with the synthesizer-based 
test signal F. The overall Lo path bandwidth is minimized by switching 
from the F + 0.495 GHz sideband to the F — 0.495 GHz sideband 
at each test set’s midband point. The proper sideband is selected by 
four stages of yia filtering. This filtering attenuates all unwanted 
tones, especially the F tone, by 140 dB with respect to the wanted 
sideband tone. This technique for the generation of the Lo signal 
eliminates the need for a second microwave source and achieves a 
coherence between the F and the F + 0.495-GHz signals that is im- 
portant in reducing phase noise. This latter consideration is discussed 
further in Section 4.4. 


3.4 Comparison circuit 


The comparison circuit is a switching network which excites the 
DUT at various test ports and. performs a rapid comparison switching 
to an internal transfer standard. In each test set, there are test ports 
with a reflection purity exceeding 40 dB used for insertion measure- 
ments and a four-configuration switching network for the measurement 


698 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1976 


of the four S-parameters of a two-port device. This latter network is 
illustrated in Fig. 4. A computer-directed calibration procedure pro- 
vides the data for the automatic removal of errors arising from. cross- 
talk and impedance deviations. Additional test ports are provided in 
the 4- to 12-GHz applique for the measurement of reflection coeffi- 
cients in three waveguide sizes. 

The rapid comparison switch allows measurements in a differential 
mode to increase accuracy. The puT is compared to the internal 
standard and the difference in loss (Maz — Ms|dB) and phase (6, — 6,) 
are recorded. Similarly, an external standard is measured, obtaining 
(Mz — Ms) and (6, — 6,). The loss and phase of the put relative to 
the external standard are derived from the simple operations: 


(Mx — Mz) = (Mx — Ms) — (Mz — Ms) dB (4) 
(0x — 62) = (6x — 6s) — (62 — Os) deg. (5) 
This technique eliminates errors from drifts in the hardware external 
to the comparison switches and from the effects of asymmetry in the 
two switch positions. The comparison rate was chosen as 30 Hz to 
achieve an effective filtering of 60 Hz ac interference. 
3.5 Interface circuitry 


Digital control of the test sets by the corms computer and computer 
monitoring of their circuit status is provided by the interface circuitry. 
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Fig. 4—Block diagram of the comparison circuit. 
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This circuitry processes all digital lines between the sets and corms, 
with manual control of each line also provided by a toggle switch. 
Computer or cotTms-generated logic signals provide control of the 
following major functions: s-x comparison switching, test-port 
configuration selection, put RF level control, and yie filter tuning 
through a p/a converter. 

The interface also monitors critical test system circuits for opera- 
tion within established limits. Detection of an out-of-limit or fault 
condition halts the measurement process and energizes a panel lamp 
indicating the condition. Interface panel meters associated with each 
of the monitored circuits are also provided as troubleshooting aids. 


IV. SYSTEM DESIGN CONSIDERATIONS 
4.1 General 


Previously, measurements of the required accuracy (see Table I) 
were made only with equipment of limited range and with painstaking 
care on a point-by-point basis. The need for measurements over ex- 
tended ranges, with measurements performed rapidly and automati- 
cally, required some special system designs. Some of the techniques 
needed to achieve frequency accuracy, high signal-to-noise ratios, 
thermal stability, and impedance purity are discussed below. 


4.2 Frequency accuracy 


Frequency accuracy is necessary to minimize the uncertainties in 
the measurement of loss, phase, and delay of networks that have large 
variations of these parameters with frequency. The delay determina- 
tion (see eq. 1) is directly sensitive to inaccuracies in frequency which 
lead to uncertainties in the quantity (FP. — F). 

Exceptional frequency accuracy is achieved in each test system by 
the generation of all signals from the synthesizer and reference signals 
from coTMs; these signals are derived from the same high-quality 
crystal source. 


4.3 Amplitude noise 


An amplitude signal-to-noise ratio exceeding 80 dB was a design 
objective to achieve a loss-measurement resolution of 0.002 dB for 
small values of loss. This was achieved through the use of a narrow 
detection bandwidth in corms (500 Hz) and through noise-suppression 
circuitry. 

The amplitude-noise suppression occurs in the level-control cir- 
cuitry which appears in the source and local-oscillator sections. These 
level circuits serve the dual function of maintaining nominal signal 
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Fig. 5—Phase-noise cancellation mechanism. 


levels and suppressing amplitude noise in a bandwidth exceeding 500 
Hz centered about the leveled carrier signals. 


4.4 Phase noise 


A design consideration of considerable importance was the ability 
to measure delay with small uncertainty. The objective was to have 
an uncertainty not exceeding 0.1 ns (rms) when measuring narrow 
filters, where the frequency aperture, /',; — Ff; ineq. (1), couldn’t exceed 
0.5 MHz. This objective is commensurate with a signal-to-phase noise 
ratio exceeding 80 dB.* 

Phase noise was reduced to accomplish this objective in two ways: 
the first by reducing the phase noise present on the microwave signal 
and the second by noise cancellation in the frequency-conversion and 
phase-detection circuitry. Phase noise is increased in the microwave 
source generation process, for it is multiplied in direct proportion to 
frequency multiplication. The amount of multiplication has been re- 
duced in each measuring system by the use of an intermediate stage 
of heterodyne conversion with a relatively phase-noise-free signal fe 
(see Fig. 2). The reduction is approximately four times or 12 dB. At 
12 GHz, where the multiplier is largest, the phase signal-to-noise 
ratio is approximately 30 dB. 

The second means of phase-noise reduction is by cancellation, of 
which there are two mechanisms. The first is demonstrated in Fig. 5. 

Phase noise in the source 6,(¢) propagates via the measurement path 
and the Lo generation path, appearing at ports 1 and 2 of the mea- 
surement-path iF mixer. The phase terms [6;; 7 = 0, 1, 2, ---, 5jare 
constants that are not important to this evaluation. The delay terms 
[r:i, 7 = 1, 2, 3] account for differences in the lengths of the various 


“A simple calculation from eq. 1 reveals the following relationship between the 
variance in 7 and the variance of the measured phase: ¢, = o9(v2)/[360(F'2 — Fi) J. 
An uncertainty, ¢, = 0.1 ns, is commensurate with os = 0.01 degree. 
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signal paths. The iF signal results from the argument differences of 
the input signals, which causes phase-noise cancellation to the follow- 
ing degree: 
On(t + 73) — On(t + 71). (6) 
The variance’ of this ir phase signal is equal to 
of, = ELOn(t + 73) — On(t + 71) = 2[Re,(0) — Bo, (7)], (7) 


where 7 = 73 — 71. When 7 = 0, there is perfect cancellation. For a 
worst-case noise distribution,’ assuming 7 S 100 ns, and an 1F band- 
width $500 Hz the pm-s/N is improved in excess of 76 dB. Phase noise 
is similarly canceled in the phase-reference-path mixer. 

The second phase-noise-cancellation mechanism arises from the use 
of a phase-reference signal. The difference in phase between the mea- 
surement signal and the phase-reference signal is detected. Phase noise 
common to both signals is rejected. This mechanism cancels phase noise 
originating in the source signal and the 0.495-GHz signals shown in 
Fig. 5. The cancellation exceeds 40 dB for the corms phase detector. 

The combined effectiveness of the above mechanisms was measured 
by deliberately adding sinusoidal phase modulation exceeding 20 
degrees peak to peak to the microwave source channel. No increase 
in the variance of a phase measurement (typically 0.01 degree) was 
detected. 


4.5 Thermal stability 


The rapid comparison technique of measurement described in Section 
3.4 allows slowly varying changes in the test system lying external to 
the comparison switches to be canceled in the difference measurements. 
Changes in the transmission parameters of components lying inside 
the switches are not canceled. The latter component count is large due 
to the need to switch among the many test ports. The primary cause 
of parameter change is temperature fluctuation. 

The temperature of the component assembly was stabilized in each 
test set by attaching it to a large aluminum plate, which in turn is 
surrounded by an insulating material. The composite structures have 
a thermal time constant exceeding 8 hours. This design has made 
thermal drift of these components a negligible source of error. The 
cables that connect the put to this circuitry are not enclosed in the 
stabilized environment and are a significant source of error (see Sec- 
tion 5.2). 


* The stochastic variable @, is assumed to be a wide sense stationary process.‘ 
All energy concentrated at the extreme edges of the 500-Hz ir bandwidth. 
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4.6 Test port impedance invariance 


By careful design, the impedance seen looking into the test ports 
used for simple insertion-loss measurements varies less than 2 percent 
from the nominal characteristic value of the coaxial line or waveguide. 
This is necessary to minimize reflection errors. With current hardware, 
it is not possible to achieve this impedance purity on the more com- 
plicated S-parameter switching networks. For measurements on these 
ports, the reflection errors are mathematically removed after the actual 
test set impedances are determined. This determination is accomplished 
through the measurement of calibration standards and a special hard- 
ware design.’ 

The S-parameter hardware shown in Fig. 4 is designed so that the 
impedances seen looking into each test port are constant, regardless 
of which parameter is measured. This invariance makes all S-parameter 
measurements referenced to the same impedances* and allows for the 
selfmeasurement of the physical impedances Z1 and Z2 seen by the 
puT. This consistency allows the transformation of the measured data 
into an S-parameter representation free of reflection errors. 

The invariance of the test port impedances in each test system is 
achieved by the use of isolators, attenuators, and the directivity of 
couplers to mask impedance changes caused by the switches. Variation 
of test port impedances of less than 0.05 percent has been achieved. 


V. MEASUREMENT ERRORS 
5.1 General 


The nature of error mechanisms must be understood to be able to 
predict uncertainties that may arise under actual test conditions. This 
understanding was obtained during the development of the two test 
systems by directly characterizing the error-producing mechanisms 
and then by observing the errors in characterizing networks of pre- 
dictable properties. Confidence was obtained from the consistency of 
the predicted errors and observed errors. The most important error- 
producing processes in the appliques are: thermal instabilities, non- 
linearities, crosstalk, second-order-noise products, and impedance- 
reflection errors. 


5.2 Thermal instabilities 


Thermal instabilities, which cause drift of the apparatus in the com- 
parison unit between calibration and device measurements, are the 


* In the S11 switch configuration, the detected loss and phase is directly related to 
S11(Z1, Z2) of the device under test; Z1 and Z2 are the test port impedances. 
Once Z1, Z2 and the four-parameter set S(Z1, Z2) are determined, it is possible 
through known transformations to obtain another representation such as S(50, 50). 
The transformed parameters are then free of reflection errors. 


AUTOMATED NETWORK ANALYZERS = 703 


largest source of uncertainty in absolute measurements. This error 
mechanism was anticipated in the system designs and measures were 
taken to stabilize the temperature of critical components. These mea- 
sures are discussed in Section 4.5. 

The paths that connect the test set connectors to the put are not 
controlled in temperature. These paths typically contain several feet 
of coaxial cable to achieve a flexible connection to the put. In the case 
of waveguide measurements, they also contain coax-waveguide trans- 
ducers.* Careful measurements have correlated almost all of the ob- 
served output data drifts to temperature changes in these external 
connecting paths. The following sensitivities in both test sets have 
been observed : 

A6/AT = 0.1°/F° per GHz (8) 
and 


AL/AT = 0.01 dB/F® per GHz. (9) 


External path temperature cycling of 1°F peak to peak has been 
observed over a 1-hour period in the air-conditioned laboratory en- 
vironment. The slowly varying thermal drifts do not cause significant 
uncertainties in the more critical differential measurements of delay 
[see eq. (1)] and loss slope. 


5.3 System nonlinearities 


The onset of nonlinearities in the test systems determines the high- 
signal-level limit of operation. Signal levels and components have been 
chosen to assure a measurement uncertainty less than 0.001 dB and 
0.01 degree due to deviations from linear detection. Devices with gain 
are measured at a reduced excitation level so that the output-signal 
level lies in the linear-detection region. 

One technique used to confirm linear performance consists of mea- 
suring a 6-dB loss unbalance between the two comparison-switch por- 
tions at the highest level and at a lower level. Adequate isolation exists 
between the comparison switch and the level-adjusting network to 
avoid errors due to impedance changes. The differences between the 
nominal 6-dB measurements are then used as a measure of errors 
caused by nonlinear behavior. The data indicates deviations from 
linear behavior are within accuracy objectives for both test sets. The 
components that limit linear behavior are the modulators in the 
heterodyne circuit. 


* The coax-waveguide transducers are reciprocal and low Joss for the four S-param- 
eter measurements. The transducers include isolators and matching attenuators for 
the simple insertion-loss measurements. 
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Fig. 6—Second-order noise product mechanism. 


5.4 Crosstalk 


Measurement uncertainties occur for puTs with high loss when 
signals can reach the detector by other paths than through the pur. 
By design, these leakage signals are more than 80 dB below the zero- 
loss-detected signal. For a put with loss Mz(dB), the uncertainty due 
to leakage or crosstalk is as follows: 


| Uncertainty (dB)| < 8.68 X 10#=s0/20 (10) 
| Uncertainty (deg) | < 57.3 X 104%2-80/20, (11) 


For example, the uncertainties for a device having 60 dB of loss could 
be 0.9 dB and 5.7 degrees. Crosstalk uncertainties are removed from 
reflection-coefficient and return-loss measurements but not from trans- 
mission measurements.* The second uncertainty term in eq. (8) ac- 
counts for this crosstalk. 


5.5 Second-order noise products 


The high degree of correlation between the rF signal and the RF 
local oscillator signal in both appliques leads to a small bias in the 
measured data. The mechanism can be understood by observing the 
simplified diagram of half of the comparison heterodyne section shown 
in Fig. 6. The harmonically generated microwave signal e; has asso- 
ciated with it a noise spectrum whose envelope is determined by the 
yi@ filters and the spectral shape of the noise phenomena. The phase 
noise spectral shape typically follows a 1/Af decay about the carrier 
frequency. The amplitude noise is typically flat and is significant only 


*The removal of crosstalk uncertainties from the transmission measurement is 
possible with additional calibration measurements and data processing. 
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in the 4- to 12-GHz test set due to the use of rwT amplifiers with noise 
figures exceeding 40 dB. 

A frequency-shifted replication of the noise spectrum of é; is present 
on the local oscillator signal e3. A noise component on é,, illustrated 
in Fig. 6 as a dashed line, interacts in the converter with a replica of 
the same component on e; to produce a coherent IF signal. The cumula- 
tive effect of all noise components produces a small signal at frequency 

= 495 MHz, leading to a bias in the measured quantities. 

A detailed analysis* of the phenomena (see Appendix) indicates that 
the average value of the ir signal amplitude must be adjusted by the 
multiplicative factor, 


[1 + Ram(r) + Rpm(r)]. (12) 


Ram(r) and Rpm(r) are the autocorrelation functions of the am and 
PM noise spectrums, respectively, that are common to the signals e; 
and e3. Similarly, the phase is adjusted by the additive factor, 


[2 Ram/pm(r) ] radians. (13) 


Ram/pm(r) is the cross-correlation function of the am and PM noise 
spectrums. The parameter (7) is the delay difference between the 
two paths to the converter shown in Fig. 4. It is therefore a function 
of the delay through the put. 

The uncertainty caused by this second-order-noise interaction is 
proportional to the changes in the above terms due to the differential 
nature of the measurements [see eqs. (4) and (5) ]. These uncertainties 
are to first order: 


Uncertainty (dB) = 8.68[Rdm(r) + Rpm(r)](rx — rz) (14) 
Uncertainty (deg) = 115[Rdm/pm(r) ](rx — 72). (15) 


The delay difference (rz — rz) is the difference in delay between the 
put and the calibrating device. The prime denotes the first derivative 
with respect to (r). 

This analysis gave the qualitative understanding of this phenomena 
that led to the use-of an intermediate stage of heterodyne conversion 
to reduce source noise (see Section 4.4). A quantitative understanding 
is prevented by a lack of knowledge of the noise-process statistics. 
Measurements made to reveal these uncertainties have shown that 
they are negligible for relatively broadband networks and for differ- 
ential measurements such as delay and loss slope. Uncertainties of 
0.01 dB and 0.03 degree were observed on the steep skirts of a 10-GHz 


* Assuming the heterodyne converter is approximated as forming the product e2-e3. 
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cavity filter with a 40-MHz bandwidth. The uncertainties decrease 
directly with frequency since the phase noise dominates (see Section 
4.4). These uncertainties are commensurate with the performance 
discussed in Section 2.3. 


5.6 Errors from impedance mismatches 


Errors are made in measuring a puT when this device, calibration 
standards, and the test set terminals have impedances that differ 
from the assumed nominal values. In the microwave test systems, 
these errors are quite significant and are fundamentally different in 
the various measurement modes. 

The simplest mode is the measurement of insertion transmission. 
The impedance mismatch uncertainty in measuring either the forward- 
or reverse-scattering parameters S2; and Si. has been discussed by 
many authors.!?)> A general relationship is: 





Siz _ 1 — poSiz — pis + pop AS 
1 — pgpi 





(16) 


S 


». 


t 
where 


M;; is the voltage-scattering parameter actually measured between 
the test set generator impedance Zg and the test set load 
impedance ZI. 

_ 49 — kg _ 41 — Kl, 
Oe Gow Re Fo SEER 


The S-parameters are those that would be obtained if the test set had 
ideal impedances Rg and Rl. 

In both microwave test systems, the insertion transmission mea- 
surements are made with relatively simple hardware. This simplicity 
allowed the impedances of the test ports to be kept within 2 percent 
of nominal. Examples of the impedance mismatch uncertainty that 


(17) 


can occur in measuring unknowns (assuming |Si2-S2i| $1; |p9| 
<0.01, |p:| S 0.01) of varying input and output return losses have 
been computed from eq. (16) and are given below: 


put Return Loss’ Max. Amplitude Error Max. Phase Error 


40 0.0035 dB 0.023° 
30 0.0072 dB 0.048° 
20 0.0191 dB 0.126° 


The impedance or reflection coefficient of a put is determined 
from one measurement of the put and one measurement on each 


* Return loss = —20log|Si;|_ 7 = 1, 2. 
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of three calibration standards. The standards are typically two 
reactive circuits (nominally an open and a short) and a resistive 
circuit (Zs) nominally matched to a transmission line characteristic 
impedance (Rs). The uncertainty that is incurred if Zs # Rs is (Ref. 2): 


Ne = 8 Pe a) (18) 
1 — pT 
Zs — Rs 

ps = Zs a Rs (19) 


M is the measured voltage-reflection coefficient and I is the reflection 
coefficient that would be measured if Zs = Rs. p, is smaller than 0.001 
for the better resistive standards. The errors arising from uncertainties 
in the reactive standards are much smaller and are ignored. Notice 
that (18) vanishes for l = +1. This occurs when the reflection mea- 
surement reduces to a differential comparison of the puT with either 
the open or short standard. Notice that the use of calibration standards 
makes the mismatch errors dependent on the quality of the calibra- 
tion standards and not upon the impedances of the test system’s 
terminals. 

Four-parameter characterizations are made with complex switching 
hardware in which the test-port impedances can vary with frequency 
up to 20 percent from the design values. The large impedance-mis- 
match uncertainties that could arise are reduced by processing mea- 
surements of the puT with measurements on calibration standards. 
The characterizations are performed in three steps. First, all S parame- 
ters defined with respect to the actual test set impedances, Z1 on port 
1 and Z2 on port 2, are determined (see Section 4.6). These parameters 
are, by definition, free of impedance-reflection errors. Next, the im- 
pedances Z1 and Z2 are measured using the same three reflection co- 
efficient standards mentioned above. Finally, these S(Z1, Z2) parame- 
ters are transformed to another representation, most commonly 
S(Rs, Rs). The reflection errors that arise in the transformed represen- 
tation are dependent upon the accuracy to which Z1 and Z2 or their 
voltage-reflection coefficients p; and p2 are characterized. Thus, the 
impedance mismatch errors are related to uncertainties in the cali- 
bration standards and not the test set impedances. 

The quantities p,; and p2 are,* from eq. (18), approximately equal 
to pei and ps2, respectively. p.1 and ps2 are the reflection coefficients of 
the resistive calibration standards R,, and R,2 used on ports 1 and 2, 
respectively. The relationships between the measured parameters 


“lof K1 t=1,2. 
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[S(Zs1, Zs2) ] and the error-free parameters [.S(Rs1, Rsz)_] are found 
by known transformations.? 


Sij(Re1, Ree) ik l= pie — Peach Pihsghs 

Si;(Ze1, Ze2) (1 a psi) (1 — ps3) Rs Roo 
=1 — pyi[ Sii(Rei, Ree) + 1) 

— pejLS5j(Rat, Rea) _ 1]. (20) 

Sit — psi + peips S53 — psjiAS 
ae psiSii = psjSij oF PsipsjAS Rai, R52 
Siis(Be1, Rez) — psilLl — Si2(Rer, Ree) ] 

+psjSij(Rsi, Reo) Sji(Rsi, Rez). (21) 





Sii(Zs1, Zs2) 





These relationships coupled with extensive experimental data were used 
to derive the uncertainty bounds given by eqs. (2) and (8). 


VI. MEASUREMENT RESULTS 
6.1 General 


Ideally, measurement accuracy should be confirmed with network 
standards selected to accentuate various error mechanisms. At micro- 
wave frequencies, no complete set of such networks exists. The mea- 
surement results given below* examine accuracy with a few network 
standards whose properties are based upon geometric dimensions. 
Additional indications of accuracy are obtained from characterizations 
of a given network under different measurement conditions that are 
expected to give the same results, and by comparing characterizations 
on both iest systems in the 3.7- to 4.2-GHz band of overlap. 


6.2 Transmission measurements 


6.2.1 Delay of an airline 


The ability to accurately measure envelope delay is of critical im- 
portance to the development of transmission networks. The delay- 
measuring uncertainties in both measuring sets were evaluated by 
measuring an airline of known electrical length, one of the few avail- 
able standards of delay. The data shown in Fig. 7 is a measurement of 
the delay of a 7-mm airline of 20.35 cm in electrical length. A fre- 
quency difference of 20 MHz was used in the on-line calculation of 
delay using eq. (1). The results show a measurement uncertainty of 
less than +0.022 ns about the nominal delay of 0.678 ns. 


* Errors caused by the failure of connectors to make good contact are not discussed. 
These errors are not significant for waveguide connections. For this reason, many 
of the networks measured for accuracy verification have waveguide connectors. 
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Fig. 7—Delay of an airline with a 20.35-cm electrical length. 


6.2.2 Insertion parameters of a waveguide filter 


The loss and delay of a narrow-bandwidth waveguide filter was 
measured on both test sets to compare their accuracy. Waveguide 
isolators were used to improve the impedances of the coax-waveguide 
transducers needed to make measurements on the 1- to 4-GHz test 
set. The same isolators were used in the measurement made on the 
4- to 12-GHz test set so that the puT would see essentially the same 
test set impedances in each measurement. The filter characteristic 
showing the passband and a portion of the stop band measured on 
the 1- to 4-GHz test set is shown in Fig. 8. Equivalent data taken on 
the 4- to 12-GHz test set (not shown) coincides point for point except 
at loss levels above 80 dB. At loss values above 80 dB, the measured 
loss is limited by crosstalk to about 88 dB. 

The difference between characterizations on the two sets in a mag- 
nified portion of the passband is shown in Fig. 9. The worst-case loss 
and delay* differences are 0.004 dB and 0.08 ns. The average absolute 
differences are 0.0012 dB and 0.019 ns. These differences are less than 
the bound on the typical uncertainty computed from eq. (3), 0.005 dB 
and 0.1 ns. 


* Delay was computed on-line using eq. (1) with a frequency difference of 1 MHz. 
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Fig. 8—Insertion loss of a waveguide filter measured on the 1- to 4-GHz measuring 
set. 


6.2.3 Measurement stability 


The insertion loss and phase of a 12-dB attenuator was measured 
at the beginning and end of a 20-minute time interval to evaluate the 
stability of insertion measurements. The attenuator was surrounded 
by thick insulation to increase its response time to temperature 
changes. Such precautions are sometimes necessary to make precise 
measurements at microwave frequencies depending on the tempera- 
ture dependence of the pur. 
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Fig. 9—Passband loss and delay of waveguide filter measurements. 


The differences in the two measurements are shown in Fig. 10. The 
loss differences are typically smaller than 0.004 dB. The phase differ- 
ences are approximately proportional to frequency as if there were a 
change in the electrical length of the attenuator or test-system cir- 
cuitry. The data indicate an equivalent change in electrical length of 
less than 0.001 inch over the 20-minute interval. This corresponds 
to a change of less than four parts in one million in the more than 20 
feet of equivalent circuit length associated with either measuring 
system (see Section 4.5). 


6.3 Return-loss measurements 


Return-loss or one-port characterizations are measured on the 4- to 
12-GHz test set on both waveguide and coaxial test ports. These 
measurements are made on the 1- to 4-GHz test set using coaxial test 
ports. Low-loss transducers can be used on either set to make mea- 
surements in other connectors. The measurement of the pur and 
calibration is under computer direction. Three calibration standards 
described by a constant magnitude and a phase that varies linearly 
with propagation coefficient are used. This description allows for the 
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correction of some perturbations from nominal values and for the 
mathematical translation of the reference plane of measurement. The 
calibration standards typically chosen are a nominal short circuit, 
open circuit, and a resistive termination. These features will be demon- 
strated below in the measurement of some precisely machined wave- 
guide networks. 


6.3.1 Inductive iris 


A symmetric inductive iris 1.7 inches wide by 0.2 inch thick was 
constructed in WR229 waveguide and measured on both test systems. 
Measurements on the 4- to 12-GHz system were made with connection 
to a high-directivity magic-tee hybrid. Measurements on the 1- to 
4-GHz system were made with connection to the internal coaxial 
hybrid through a lossless waveguide to coax transducer. Calibration 
was performed with three waveguide standards, a short, an offset 
short, and a match. A different set of standards was used on each test 
system. Figure 11 gives the measurement results on the 1- to 4-GHz 
system. Data obtained from the 4- to 12-GHz system (not shown) 
exhibit a mean absolute difference of 0.01 dB and 0.08 degree from 
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Fig. 10—Change in insertion loss and phase of a 12-dB attenuator after a 20- 
minute time interval. 
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Fig. 11—Insertion loss and phase of a symmetrical inductive window 1.700 inches 
wide by 0.2 inch thick. 


the 1- to 4-GHz system data, with peak differences of 0.03 dB and 
0.15 degree. The solid data points are derived from theoretical com- 
puter computations based upon the mode-matching technique.* These 
differences are less than the uncertainty of 0.037 dB and 0.24 degree 
computed from eq. (2). Experimental results in the next section re- 
veal that a major contributor to the uncertainties in these measure- 
ments is the deviations of the match from a zero reflection coefficient. 


6.3.2 Offset short 


A short circuit displaced 0.4090 cm in WR90 waveguide was mea- 
sured as a circuit standard. This circuit has a reflection coefficient of 
unit magnitude (0 dB) and an angle of —90 degrees in the middle 
of the 10.7- to 11.7-GHz microwave radio band. The calibration 
standards used for the measurement were a zero length short, a match, 
and a short offset by 0.8197 cm. The latter calibration standard has a 
reflection coefficient of unit magnitude and an angle of 180 degrees in 
midband. The match has a return loss exceeding 60 dB. The measure- 
ment of this circuit should produce the worst-case uncertainties since 
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Fig. 12—Return loss of an offset short in WR90 waveguide. 


its reflection coefficient has the greatest distance on the Smith chart 
from the reflection coefficients of the calibrating standards. 

The measurement results are shown in Fig. 12. The measurement 
reference plane was mathematically moved to the plane of the circuit’s 
short circuit by an on-line program. The nominal characterization at 
this displaced plane is a reflection coefficient of unit magnitude and 
constant 180-degree phase angle. Two characterizations are presented, 
one with the tapered absorber in the match in its nominal position 
and a second with the absorber moved 0.83 em (~2/4) within its 
electroformed waveguide. The changes in the characterizations reveal 
that deviations of match from a zero reflection coefficient are a signifi- 
cant source of measurement uncertainty. The worst-case deviations 
are less than the uncertainties of 0.052 dB and 0.34 degree computed 
from eq. (2). 


6.4 Four-parameter characterizations 


Four-parameter or two-port characterizations are measured on both 
microwave test sets using automatically switched hardware and a 
computer-directed calibration technique. These features were briefly 
discussed in Sections 2.2 and 5.6. Voltage-scattering parameters (S 
parameters) are computed in real time. The calibration-standard 
models used assume a constant magnitude and a phase that varies 
linearly with propagation coefficient. This description allows for the 


AUTOMATED NETWORK ANALYZERS 715 


0.01 


© 
oO 
nN 


LOSS IN DEGREES 
° 
ro) 
@ 





0.04 
1.0 2.0 3.0 4.0 


FREQUENCY IN GHz 


Fig. 183—S21 phase of a 50-ohm airline. 


correction of some perturbations from nominal values and for the 
translation of the reference planes of measurement. Some of these 
features will be demonstrated in the measurements below. Parameter 
representations other than S parameters are available through process- 
ing following the measurements. 


6.4.1 Airline characterizations 


An airline is a two-port network whose characteristics are precisely 
known. Errors can arise in characterizing an airline if the test set 
impedances deviate significantly from nominal. Impedance deviations 


0.05 


0.10 


0.15 


PHASE IN DEGREES 


0.20 





0.25 


1.0 2.0 3.0 4.0 
FREQUENCY IN GHz 


Fig. 14—S21 magnitude of a 50-ohm airline. 
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Fig. 15— S11 magnitude of a 50-ohm airline. 


in excess of 20 percent are observed in the test sets. These errors are 
avoided by the calibration method discussed in Section 5.6. 

A 30-cm-long 14-mm airline was characterized on the 1- to 4-GHz 
measuring system to illustrate these points and to demonstrate the 
accuracy obtained. The S21 measurements are plotted in Figs. 13 and 
14 on a square-root-of-frequency scale. The reference planes of mea- 
surement were translated to the center of the airline so that the phase 
data is the difference between the measured airline and an ideal airline. 
Observe that the data indicates characteristics nearly linear* with re- 
spect to the square root of frequency as expected from theoretical 
models, taking skin losses into account.? The maximum deviation from 
linear is 0.004 dB and 0.05 degree. The typical deviation derived from 
eq. (3) is predicted to be less than 0.005 dB and 0.04 degree. An S21 
measurement was also made without calibration to demonstrate the 
accuracy Improvement provided by computationally removing the test 
set impedance-reflection errors. Deviations from a square-root-of-fre- 
quency characteristic in excess of 0.1 dB and 0.7 degree were observed. 

The magnitude of S11 (dB) is shown in Fig. 15. The uncertainty of 
this measurement is dominated by deviation of the resistive calibration 


* The solid line in Fig. 14 was extrapolated from resonant Q measurements made 
at 250 MHz. 
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standard from 50 ohms, corresponding to a reflection coefficient un- 
certainty of 0.005 (46 dB). 


6.4.2 Transistor characterization 


The 1- to 4-GHz measuring system is frequently used for the charac- 
terization of microwave transistors to be used in amplifier circuits. 
Their small size and impedances that vary considerably from 50 ohms 
make them particularly difficult to characterize accurately. Two 
characterization techniques have been used. The first technique uses 
simplified handling by calibrating with large coaxial standards at the 
interface of the transistor jig and test set. The defined reference planes 
are then translated to the jig-transistor interface making the approxi- 
mation that the jig is an ideal 50-ohm transmission line. The second 
technique uses calibration standards of the same size scale as the 
transistor inserted at the jig-transistor interface. The parasitic circuit 
elements of these standards are significantly smaller than those as- 
sociated with the transistor Jig. A detailed characterization effort is 
under way so that a definitive accuracy evaluation can be made. 

Measurements of S11 using the two techniques to characterize a 
beam-leaded transistor mounted on a 0.1-inch-square substrate are 
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Fig. 16—Input reflection coefficient of a beam-leaded microwave transistor. 
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Fig. 17—Forward scattering parameter of a beam-leaded microwave transistor. 


shown in Fig. 16. The top curve is the characterization using the 
transistor-sized standards. This curve is free of the impedance cycling 
with frequency that is evident in the second curve measured with the 
coaxial standards. This cycling cannot be attributed to the transistor 
and is found in measurements on the jig with the transistor replaced 
with a through connection. 

The measurements of S21 using the two calibration techniques is 
shown in Fig. 17. The measurements using the coaxial standards cycle 
about the measurements using the transistor-sized standards, ending 
up 0.5 dB lower at 4 GHz. These results reveal that measurement 
accuracy can be improved by calibration directly at the interface of the 
put and the test set. 


Vil. SUMMARY 


The development of two general-purpose network analyzers covering 
the frequency range of 0.9 to 12.4 GHz has been discussed. These net- 
work analyzers were developed to provide flexible, fast, and accurate 
circuit characterizations, which are necessary feedback in the circuit 
design process. 

Automatic control and highly accurate baseband detection is pro- 
vided by the Computer Operated Transmission Measuring Set. The 
techniques used to extend baseband-measurement capabilities to micro- 
wave frequencies while retaining accuracy were discussed. The most 
prominent sources of measurement errors were quantitatively ex- 
amined. Measurements results that confirm accuracy were presented. 
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APPENDIX 
Analysis of Second Order Noise Products 


Rice’s representation (Ref. 4, p. 374) of a stochastic process is con- 
venient for the analysis of the heterodyne process modeled by Fig. 6. 


€; = sin (w,t) + n(t) sin (w.t) + m(t) cos (w.6). (22) 


The first term is the normalized carrier signal. The second term is in 
phase noise which can be thought of as am noise for small n(t). The 
last term is in quadrature noise which can similarly be thought of as 
phase noise. The stochastic variables are assumed to be stationary 
and to have the following properties.* 


E{n@] = Elm(t)] = 0. (23) 
E{n(t + 7)n(t)] = Ram(r) = Ram(—7). (24) 
Elm(t + r)m(t)] = Rpm(r) = Rpm(—7). (25) 


Elm(t + 7)n()] = Ram/pm(r) = —Ram/pm(—7). (26) 


The delayed signal e2 and the up-converted signal e; are given by 


€2 = sin (wt + wer) + n(t + 7) sin (wet + w-7) 
+m(t + 7) cos (wet + wer). (27) 
é3 = 2sin (w. + A)ét + 2n(Z) sin (w, + A)t 


+2m(t) cos (we + A)t. (28) 
The output converter is approximated as forming the product e4= e2: és. 


és = cos (At — w,7) 
+[n(t + r)n(t) + m(t + 7)m(t)] cos (At — w-7) 
+[m(t + r)n(t) — n(t + 7)m(t)] sin (At — 7) 
+[n(t + 7) + n(t)] cos (At — w-7) 
+[m(t + r) — m(t)] sin (At — 7). (29) 


Filterable tones falling at 2w, have been omitted from eq. (29). 
Taking the expected value of e, and substituting the properties of 
eqs. (23) through (26) yields, 


E(es) = cos (At — wr) + [Ram(r) + Rpm(r)] cos (At — wer) 
+2Ram/pm(r) sin (At — wet). (80) 


E(es) = [1 + Ram(r) + Rpm(r)] cos [At — wer — 2Ram/pm(r) J. 
(31) 


* E[-] denotes expected value. If Ram(r) = Rpm(r), then the noise process also 
models an additive as well as a modulating stationary process.‘ 
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The amplitude and phase of the mean rF signal e, are seen in eq. (31) 
to be modified by the noise present on the microwave signals. 
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An automatic speaker-verification system accessed by test customers 
from their own telephones over dialed-up lines has been evaluated. The 
test population consisted of over 100 male and female speakers who called 
up nominally once each working day over a period of five months. The 
operation of the system is based on a set of functions of tume obtained from 
acoustic analysis of a fixed, sentence-long utterance. These functions are 
compared with stored reference functions to determine whether to accept 
or reject an identity claim. The system is implemented on a NOVA 
800 laboratory computer. Telephone line access to the computer is via a 
data set hookup. Identity claims are made by keying an identification 
number on a Touch-Tone® dial. Instructions and responses to the 
customer are made by means of a programmed voice-response system. 
Reference data was computed off-line and updated with the analysis data 
of accepted utterances. The evaluation indicated an error rate of approxi- 
mately 10 percent for new customers and approximately 5 percent for 
adapted customers. 


Il. INTRODUCTION 


Speaker verification is the authentication of an individual’s claimed 
identity by analysis of his spoken utterances. Research on an auto- 
matic system for speaker verification at Bell Laboratories has been 
reported in previous papers.'~* The system is based on an acoustic 
analysis of a fixed, sentence-long utterance resulting in a function of 
time or contour for each feature analyzed. Features selected for anal- 
ysis in previous evaluations have included pitch, intensity, the first 
three formants, and selected predictor coefficients. The system com- 
pares the set of sample contours obtained from an unknown individual 
with the set of reference contours corresponding to the identity claimed 
by that individual. If the comparison results in an overall measure of 
dissimilarity which is smaller than a predetermined threshold, the 
identity claim is accepted. Otherwise, it is rejected. 
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Previous evaluations of the system have concentrated on investi- 
gating features to be analyzed and developing comparison procedures 
to make the system as effective as possible in terms of reducing overall 
error rate. The speech samples used in these evaluations were collected 
by wideband recording of male speakers in a sound booth. The record- 
ings were selected and edited to eliminate botched utterances and non- 
speech acoustic events. From the outset, however, the intent has been 
to provide a completely automatic system which could operate via 
dialed-up lines from telephones on the user’s own premises and to 
include both male and female speakers in the user population. The 
purpose of the evaluation described in this paper was to determine 
how well the system would operate under these broadened, ‘‘real- 
world’”’ conditions. 

There are several ‘real-world’ difficulties which are expected to 
be adverse to system performance. First, there are the uncontrolled 
and degraded environmental and transmission conditions encountered 
during the recording of sample utterances. Environmental conditions 
involve acoustic background noise and disturbances generated at the 
user’s end and by telephone equipment. Transmission conditions in- 
volve signal modification over dialed-up telephone lines. Telephone 
transmission is nominally over a band from 300 to 3000 Hz. The roll- 
off at 300 Hz may be gradual due to the attenuation characteristics 
of the carbon-button transmitter, or quite sharp due to the attenuation 
characteristics of repeaters in some toll lines. Moreover spectral and 
phase distortions and variations are likely to be encountered. 

The second class of problems is largely behavioral. For example, 
can a stable and adequate initial reference file be established based on 
a small number of sample utterances collected in one sitting? Also, 
can both day-to-day and long-term variations in speaking behavior 
be tolerated and tracked, and can reference files be updated to reflect 
these changes in behavior? 

Since the principal goal of this evaluation was to study the effect 
of these “‘real-world’”’ conditions rather than to achieve optimum per- 
formance, the system was made more tractable by using only pitch 
and intensity features for analysis. 


ll. SYSTEM OPERATION 


Although the operation of the system has been described in previous 
papers, it will be outlined again here, with departures pointed out, to 
provide a basis for discussion of the present evaluation. 

Figure 1 provides an outline of the system operation in the form of 
a block diagram. The entire system has been implemented in soft- 
ware on a Data General NOVA 800 laboratory computer. The two 
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Fig. 1—System-operation block diagram. 





inputs to the system are an identity claim and a sample utterance. 
The sample utterance used for purposes of evaluation is the all-voiced 
sentence ‘‘We were away a year ago.’”’ A marked interval is provided 
for input of the utterance. The input is subjected to 900-Hz, low-pass, 
analog filtering by two sections of a Rockland 1520 Dual Filter with a 
combined roll-off of 48 dB/octave. The filtered input is digitized at 10 
kHz by means of a 12-bit analog-to-digital converter and stored on 
disk. The digitized input is then scanned forward from the beginning 
of the recording interval and backward from the end to determine 
the beginning and end of the actual sample utterance. The end-point 
detection is accomplished by means of an energy calculation. The de- 
limited portion is subjected to feature analysis which in this imple- 
mentation consists of a pitch-and-intensity analysis. Pitch analysis is 
accomplished by means of the time-domain parallel-processing tech- 
nique of Gold and Rabiner,’ modified by Rabiner for application to 
telephone speech and extension to female talkers. A pitch period value 
is obtained every 10 ms through the course of the utterance with 
resolution to 100 us, the sampling period. The resulting pitch contour 
is smoothed nonlinearly to bridge across unvoiced gaps and to diminish 
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Fig. 2—Smoothed pitch period and intensity contours obtained from analysis of 
a sample of the test utterance. 


the effect of singular values.* The contour is further subject to 16-Hz 
low-pass smoothing. In addition, an intensity or energy measurement 
is calculated every 10 ms through the course of the utterance to obtain 
an intensity contour. This contour is also subject to 16-Hz low-pass 
filtering and, in addition, is normalized to the peak intensity measure- 
ment resulting in a contour of relative intensity values. A typical set 
of pitch-and-intensity contours is shown in Fig. 2. 

These contours comprise the basic patterns for verification. They 
are compared with a corresponding set of reference patterns associated 
with the claimed identity. The reference patterns are obtained by 
averaging and combining a set of patterns obtained from sample 
utterances of the person whose identity is claimed. (The referencing 
process has been described in Ref. 2.) Before comparing the sample 


“The pitch detector modifications, especially the nonlinear smoothing, were 
largely motivated by a preliminary study of the effects of telephone transmission on 
automatic speaker verification by O. M. M. Mitchell.§ 
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and reference contours, an additional operation, time registration, is 
carried out. In this operation, the events of the sample contour are 
brought into the best possible registration with corresponding events 
of the reference contour by replotting the sample contour versus a 
modified function of time. This step is necessary to account for the 
normal, expected variations in speaking behavior observed in the 
repetition of a sample utterance by the same speaker. In previous 
implementations, this operation was accomplished by means of a 
method of steepest ascent.! In the present implementation, a dynamic 
programming technique is used. The technique is similar to those 
described by Sakoe and Chiba, Itakura, and Ellis.*-8’ The intensity 
contour is the guide contour for the procedure. The sample intensity 
contour is linearly stretched or compressed to the normalized length 
of the reference intensity contour. Then a distance is calculated be- 
tween the ith point (or set of points) in the sample contour and the 
jth point (or set of points) in the reference contour for each 7 and j. 
The dynamic programming algorithm is used to find the path of least 
accumulated distances through the matrix of distances {d;}. The 
optimal path ¢ = I(j) 7 = 1, ---, N determines the warping function 
required to replot the sample contour registered to the reference con- 
tour. A number of constraints are imposed so that the resulting path 
does not deviate excessively from the path of no warping 7 = j. The 
warping function obtained for the intensity contour is also applied 
to the pitch contour. Time registration of the intensity contour is 
illustrated in Fig. 3. 

Following registration, the pitch and intensity contours are divided 
into 20 equal-length segments, as shown for intensity in the bottom 
panel of Fig. 3. In each segment, a set of measurements is applied to 
both the sample and reference contours and a squared difference is 
calculated specifying the dissimilarity between these contour segments 
for each measurement. The squared difference for each measurement 
and segment is weighted inversely by a variance which is calculated 
from the set of sample contours used to construct the reference (see 
Section 2.2). The effect of the variances is to weight most heavily 
those segments in which a particular measurement was most consis- 
tent over the set of sample contours comprising the reference. A 
distance for each measurement is calculated by summing the weighted 
squared differences over the 20 segments of a contour. In addition to 
four distances for each contour, based on segment-by-segment mea- 
surements, there is also a distance based on the overall cross-correlation 
of sample and reference contours. There are also a distance based on 
the cross-correlation of the pitch and intensity contours and two dis- 
tances based on the amount of warping required to register the sample 
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Fig. 3—Time registration of sample and reference intensity contours. The reference 
contour is plotted using a dotted line in both the top and bottom panel. The sample 
contour, solid line, is shown with its end points aligned to the reference contour, 
before internal registration in the top panel, and after internal registration in the 
bottom panel. 


contours to the reference contours. The overall distance is obtained 
by a simple average over the entire set of 13 individual distances or 
the average over a subset of these distances selected a priort for each 
speaker. The speaker-dependent distance-selection technique is de- 
scribed in Ref. 3. Finally, the overall distance is compared with a pre- 
determined threshold to determine whether to accept or reject the 
identity claim. 


2.1 Experimental setup and typical transaction 


A block diagram of the experimental setup is shown in Fig. 4. The 
basic elements are the customer’s phone with a Touch-Tone dial, a 
data set (Western Electric 407A), an analog-to-digital converter, and 
a Data General NOVA 800 computer, in which reside both the auto- 
matic speaker-verification programs and a programmed voice-response 
facility which is used to provide instructions and responses to the 

‘customer during each transaction. The voice-response system, which 
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Fig. 4—Experimental setup block diagram showing the principal components of apparatus linking the customer’s telephone and the 
computer. 


uses ADPcM-coded speech message units, has been described in a paper 
by Rosenthal et al.° 

A typical transaction begins with a customer dialing up the data set 
from his own premises via the local telephone exchange. The data set 
acknowledges the incoming call and provides an interrupt to the cru, 
which initializes an identifying voice-response message and the in- 
struction to the customer to dial in his identification number via the 
Touch-Tone dial on his phone. The decoded identification number 
constitutes the customer’s identity claim. The customer is then re- 
quested to speak the test utterance within a tone-marked recording 
interval and the input utterance is then processed as described in the 
previous section. (There are several conditions that can occur during 
the end-point detection or pitch-and-intensity analysis indicating a 
bad recording that will cause the customer to be instructed to repeat 
the test utterance.) Following the analysis and comparison, the cus- 
tomer is advised of the decision to accept or reject his identity claim. 
A rejection, of course, constitutes a system error, a false alarm. The 
entire transaction from recording to verification response takes 20 to 
30 seconds. A large fraction of this time is consumed by the software 
pitch detector. A breakdown of the computation times is shown in 
Table I. The major decision parameters for each transaction are ap- 
pended to a log file set up for each customer. 

One hundred four “customers,’’ approximately evenly divided be- 
tween adult males and females, participated actively in the evaluation 
by calling in nominally once each working day over a period of five 
months. These were all native American or Canadian speakers of 
English with no overt speech defects. The customers were instructed 
to speak the test utterance naturally and consistently from day to day. 


2.2 Reference information 


The establishment and updating of reference information is an 
important element of the system. On the customer’s initial call to the 
system he is requested to provide five recordings of the test utterance 
with approximately 10 seconds between each recording. These utter- 
ances are analyzed and the analysis data are stored in disk files. These 
data are used to construct the customer’s initial reference file. The 
actual reference construction is carried out off-line during nonoperating 
hours and requires about 10 minutes per customer. 

The following operations are included in the construction of a refer- 
ence file: 


(i) Reference contours. A set of reference pitch and intensity con- 
tours is constructed. The intensity contours obtained from up 
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Table |— Processing times for a single transaction 


: Time 
Operation (s) 
Intensity computation & end-point detection 4 
Pitch analysis 5.5 (Typical) 
(3.5 s/s) 
Time registration 3.5 
Comparison & decision 2 
15 
System overhead 10(Approximate) 
Total 25 


(2) 


(222) 


(wv) 


For 


to 10 customer sample files are mutually registered and aver- 
aged to provide a reference intensity contour. The correspond- 
ing sample pitch contours are mutually registered and averaged 
using the same warping parameters obtained for the reference 
intensity contour. 
Distance weights. For each kind of measurement made on the 
contours, a variance is calculated over the sample contours 
used to construct the reference contours in the form 
Bea Sa cen 

TNS, Sin T5)"; . 
where s;, is the jth measurement on the nth sample contour, 
r; 1s the jth measurement on the reference contour, and N is 
the number of sample contours. 
Measurement selection. The subset of the original measurement 
set which is most effective in separating the overall distance 
distribution of customer and impostor sample utterances is 
found. Each customer sample file and a sample file from each 
of 30 different customers of the same sex are used to provide 
the customer and impostor distributions. 
Threshold computation. Estimates of equal-error thresholds for 
both the overall distance based on all measurements and the 
overall distance based on selected measurements are computed. 
The same sample customer and impostor files used in measure- 
ment selection are used to estimate the thresholds. 


the initial reference file, the operating threshold is set at 1.5 


times the estimated equal-error threshold. This is done to compensate 
for the fact that the sample files used to estimate the threshold coincide 
exactly with the sample files used to construct the reference contours. 
Moreover, these sample files are obtained from utterances collected 
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in one session. For these reasons, the measurements and distances 
extracted from this set of sample files are expected to be highly cor- 
related and are not likely to adequately reflect the expected variation 
over a series of independent trials. The factor of 1.5 for augmenting 
the threshold is somewhat arbitrary and was arrived at by inspection 
of the data from a preliminary experiment. The accept-reject criterion 
following the initial reference is based on all measurements rather than 
a selected set of measurements. This, again, is done because the initial 
sample data files do not adequately reflect the expected range of vari- 
ation to provide a stable selection of measurements. 

Following establishment of the initial reference, the customer calls 
in nominally once a day. The data files for each trial in which the 
customer’s claim is accepted are saved. When five of these files are 
accumulated, the reference file is updated. The reference file is up- 
dated a second time when five additional ‘‘accepted’”’ customer sample 
files are accumulated. For the fourth reference and thereafter, 10 
additional customer files must accumulate. For the fourth and succes- 
sive references, there are a total of 25 customer files available of which 
10 are used to construct the contours and all are used to select mea- 
surements and calculate the thresholds. Following reference construc- 
tion, the oldest 10 customer files are deleted, so that the maximum 
number of sample files per customer allowed to accumulate is 25. Also, 
from the fourth reference and thereafter, the operational accept/reject 
criterion switches over to selected measurements and the threshold 
is allowed to adapt. The adaptation mechanism is as follows: at each 
trial for which the customer claim is accepted, if the overall distance 
is greater than 75 percent of the current threshold, the threshold is 
increased by 10 percent; if the overall distance is less than 40 percent 
of the current threshold, the threshold is decreased by 10 percent. 
The reference updating procedures are summarized in Table II. 


Table I! — Reference updating procedures 


No. of Utterances Analyzed 


ae No. Used Threshold Setting 
Total to Construct 
Ref. Contours 
1 5 5 1.5 X EET (all measurements) 
2 10 10 1.5 X EET” (all measurements) 
3 15 10 1.3 X EET* (all measurements) 
4+ 25 10 EET* (measurements selected and allowed 


to adapt to the next update) 


*psET = Equal-error threshold. 
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2.3 Results 


The results are described in terms of the two types of error that 
can occur: rejecting a (legitimate) customer claim and accepting an 
(improper) impostor claim. The customer-rejection rate was calcu- 
lated simply from a tabulation of the number of rejections experienced 
by the customer population over the five-month-long period of evalua- 
tion. For each trial, overall distances were calculated both for all 
measurements and for selected measurements. However, the opera- 
tional accept/reject criterion is based on all measurements for the 
first three references and on selected measurements thereafter, as 
described in Section 2.2. Excluded from this tabulation were trials in 
which it was known that an unauthorized person used an identification 
number, trials in which the customer deliberately altered his utterance, 
trials during which the analog-to-digital recording system was operat- 
ing defectively, and trials in which a faulty end-point-detection pro- 
gram produced misaligned analyses. (The latter two situations were 
quickly rectified.) Also excluded were the trials from customers whose 
first references failed. A first-reference failure is defined as one in which 
a new customer obtained three rejections before obtaining the five 
acceptances necessary for the first update. Approximately 20 or 25 
percent of the initial references failed in this way. When this occurred, 
the customer was requested to provide a new sample of five utterances 
with which to construct the initial reference. Invariably this new 
initial reference posed no problems. 

A typical customer history is plotted in Fig. 5. Each point represents 
the overall distance for a particular trial. The current threshold is 
plotted as a broken horizontal line. An error occurs for each trial in 
which a point lies above this line, that is, where the overall distance 
exceeds the threshold. Three such errors occurred over the 76-trial 
history of this customer. The reference update in effect for each series 
of trials is indicated by the numbers above the horizontal axis. The 
threshold was allowed to adapt trial-by-trial following the fifth refer- 
ence. The general trend for both distance and threshold is an initial 
elevation followed by a levelling off after 20 or 25 trials. This general 
behavior is expected and relatively easy to track. The occasional large 
trial-to-trial variations are not easy to track even with adaptation. 

The overall results of the evaluation expressed in terms of average 
error rates are shown in Table III. Three columns of figures are shown 
under “reject customer.”’ The ‘‘operational”’ criterion, as already men- 
tioned, uses all measurements for the first three references and selected 
measurements thereafter. With the operational criterion, the average 
rejection rate over all customers and approximately 4500 trials is 
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Fig. 5—Individual customer history showing the trial-by-trial overall distance 
throughout the period of the evaluation. Also shown is the operational threshold and 
reference updates through the same period. 


approximately 9 percent. The median is somewhat less, approximately 
7 percent. The error rates for female customers are consistently slightly 
higher than those for male customers. However, a statistical hypothesis 
test, derived by a likelihood ratio based on normal population dis- 
tributions, indicated that the hypothesis of identical populations can- 
not be rejected at any reasonable level of significance. Under the 
assumption of identical populations, observed differences in error rates 
could be expected with a frequency as great as 20 percent. Note that 
the operational error rate is consistently lower than either the error 
rate for all measurements or for selected measurements. This is be- 
cause the all-measurement criterion is preferable for early references 
and the selected-measurement criterion is preferable for later ones, 
as shown in Fig. 6. The left half of this figure shows the average reject- 


Table IIl— Average error rates 
Reject Customer Accept Impostor 
Customers Opera- All Selected | Opera- All Selected 
tional Measure- | Measure- tional Measure- | Measure- 
ments ments ments ments 
es 8.20 9.46 10.65 8.78 8.91 6.06 
ae 9.92 10.37 12.2 10.92 10.94 7.54 
48 
Total 8.99 9.88 11.33 9.72 9.79 6.70 
(104) 
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Fig. 6—Average error rates vs reference update. Reject-customer rates are plotted 
on the left and accept-impostor rates on the right. The thick line in each case shows 
the error rate under the operational criterion. 


customer rate as a function of reference update. The rates for the all- 
measurement criterion and the selected-measurement criterion are 
plotted separately and the combined operational strategy is indicated 
by the thick line. The general trend for both criteria is the same as 
that observed for the history of a single customer shown in Fig. 5: 
a high initial error rate followed by a levelling off. It can also be seen 
that the operational strategy produces the minimum error rate. Over- 
all, at early reference update stages, the reject-customer rate is of the 
order of 10 percent, while in later stages, the rate approaches 4 percent. 

The second part of the error analysis, the tabulation of accept- 
impostor data, is accomplished differently. It is not practical to pro- 
vide access to the system to a separate population of impostors for 
this purpose. Instead, a procedure was set up to systematically cross- 
compare selected sample utterances of each customer with the refer- 
ences of every other customer of the same sex. The sample file for 
every tenth accepted utterance of each customer is selected for this 
purpose, and the results tabulated. This is accomplished during the 
same off-line period used to update the reference files. The error rates 
shown in the right half of Table III are based on approximately 16,000 
comparisons. The accept-impostor rates are approximately 10 percent 
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for the all-measurement criterion and 7 percent for selected measure- 
ments, indicating a clear advantage for the selected-measurement 
criterion. Unfortunately, the operational criterion produces no such 
improvement. The reason for this is seen in the right half of Fig. 6, 
where the accept-impostor rate is plotted versus reference number. 
Most of the advantage for selected measurements occurs for customer 
references in the early stages of updating. At these early stages, how- 
ever, the operational criterion is the all-measurement one. In the latter 
stages of stable references there is no clear advantage for either all 
measurements or selected measurements. In these latter stages, as 
was the case for the reject-customer rate, the accept-impostor rate is 
approximately 4 percent. 

In summary, error rates of the order of 4 or 5 percent, for both 
reject-customer and accept-impostor, are obtained for customer refer- 
ences in advanced stages of reference updating. The overall error 
rates are about twice as high because of the adverse effect of errors 
occurring at early stages of reference updating. 

Another interesting question is the effect on the error rates of vary- 
ing the threshold. To determine the effect of threshold variation, all 
the customer log files were scanned by varying the actual thresholds 
in steps of 10 percent and tabulating the number of errors at each 
step by comparing the actual overall distances with the varied thresh- 
old. The results are shown in Fig. 7. All the actual thresholds are 
normalized to 1. Thresholds are plotted for the operational accept/ 
reject criterion. At the normalized threshold value of 1, the error rate 
is approximately 9 percent, the same value shown in Table III. The 
plot of accept-impostor error rate as a function of threshold was ob- 
tained by comparing selected customer samples left on file at the 
end of the experiment with customer references. (Since the comparison 
data is from the period at the end of the experiment, most references 
were in an advanced updating stage and the error rate at the normalized 
threshold value of 1 is approximately 5.5 percent, considerably less 
than the 9.7 percent rate obtained throughout the entire period of the 
experiment, as shown in Table III.) It is possible to get a feeling for 
the amount of tradeoff obtainable when the threshold is varied. For 
example, if the threshold is set for a reject-customer rate of 4 percent, 
the corresponding accept-impostor rate is approximately 15 percent. 
Conversely, if the threshold is set for an accept-impostor rate of 4 
percent, the corresponding reject-customer rate is approximately 11 
percent. 

Finally, it is of interest to survey individual error rates to get a 
feeling for the range of performance over the customer population. 
Histograms of individual error rates have been plotted in Fig. 8. The 
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Fig. 7—Error rate vs normalized threshold. The effects are given of threshold 
variation on the reject-customer and accept-impostor rates using the operational 
accept/reject criterion. 


top half shows the distribution of individual reject-customer rates 
under the operational criterion. About 80 percent of the customers 
have error rates less than 15 percent. There is, however, a small frac- 
tion of the population with excessively large rates of rejection. Some 
of these customers have been identified as special cases and will be 
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Fig. 8—Histograms showing the distribution of individual error rates over the 
104—customer population. The top half shows the reject-customer distribution while 
the bottom half shows the accept-impostor distribution. 


discussed in the next section. In the bottom half of the figure, a 
similar histogram is shown for the accept-impostor rates under the 
selected distance criterion. This distribution is somewhat tighter, with 
fewer outliers than the reject-customer distribution. 


IH. DISCUSSION 


As stated in Section I, the purpose of this evaluation was to deter- 
mine how well the system could perform under “real-world” condi- 
tions. First, ‘real-world’? conditions make it difficult for reference 
files to adequately represent and keep pace with normal variations in 
speaking behavior. To avoid encumbering the customer, the initial 
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reference file is constructed from a small number of highly correlated 
utterances collected in one sitting. Two questions are (7) how adequate 
is this initial reference file, and (72) can it be satisfactorily updated and 
adapted to track both trial-to-trial and long-term variations in speak- 
ing behavior? In contrast, in previous laboratory evaluations, refer- 
ence files were constructed from independent samples generally span- 
ning a relatively long period of time with a distinct set of test samples 
spanning the same period of time. 

Second, the recordings are obtained over dialed-up lines from the 
customer’s own telephone and are therefore uncontrolled and de- 
graded in comparison with the carefully executed wideband record- 
ings of previous laboratory evaluations. 

In a previous laboratory evaluation, the equal-error rate using a 
pitch-and-intensity analysis was approximately 6 percent with the 
all-measurement criterion and 3 percent with the selected-measure- 
ment criterion.* An error rate of approximately 5 percent would be 
considered quite satisfactory for the “real-world” evaluation. In fact, 
an overall error rate of approximately 9 percent was obtained. As antici- 
pated, the error rate varied considerably over the course of the cus- 
tomers’ access and reference-update history. For well-established and 
adapted customer references, an error rate of approximately 4 percent 
was obtained, which is quite acceptable. However, the error rate 
of 15 or 20 percent obtained for initial customer references is 
unsatisfactory. 

This “start-up” or initialization problem is in fact compound. In 
the first place, the number of sample utterances available during the 
first few reference updates is only marginally adequate to calculate 
reliable reference data. This is especially true for the initial reference in 
which the same five utterances are used to construct contours as well 
as to calculate weights and thresholds. It is preferable to have an in- 
dependent set of utterances with which to calculate weights and thresh- 
olds. The second part of the “start-up” problem is the customer’s 
talking behavior. More than likely there will be large variations in 
talking behavior from trial to trial through the ‘‘warm-up”’ period of 
early trials until a stable or habituated talking pattern is éstablished. 

As just mentioned, the initial reference is a special case because the 
five sample utterances used to calculate the reference data are collected 
in one session. Generally, the talking behavior from utterance to 
utterance in this session will be highly correlated. The reference data 
calculated from this set of utterances will therefore be quite “tight,” 


* With a predictor coefficient analysis added to the pitch-and-intensity analysis, 
this same evaluation yielded a 3-percent and a 1-percent error rate for all and selected 
measurements, respectively. 


SPEAKER VERIFICATION 739 


encompassing only the limited range of behavior of that session. The 
utterances in subsequent sessions, however, can be expected to vary 
considerably from those in the initial session and from each other 
due to variations In behavior during the ‘“‘warm-up” period. Without 
compensation, then, a large customer-reject rate can be expected 
for the trials immediately following the initial reference. The compensa- 
tion that was attempted was to augment the calculated threshold by 
50 percent. However, it is clear that this was not sufficient since 20 
or 25 percent of the initial references “failed,’’ as described in the 
previous section, and the elevated rejection rate was 14 percent or 
15 percent at this stage, as shown in Fig. 6. Additional augmentation 
could well be tolerated since the accept-impostor rate at this stage is 
only 6 percent. Thus, the threshold could have been adjusted to yield 
approximately a 10-percent equal-error rate. 

The situation is quite different at the first reference update, reference 
no. 2. At this stage, the sample utterances used to calculate the refer- 
ence data consist of the five utterances from the initial session plus 
five more obtained in succeeding trials. The latter, as just mentioned, 
can be expected to vary considerably from the initial set and from 
each other. The result is a set of reference data which is likely to be 
considerably “looser” than the initial set. Consequently, there is a 
considerably reduced customer-reject rate together with a consider- 
ably increased impostor-accept rate. The same situation still holds 
true at the second reference update when an additional five utterances 
are added to the reference data. In both of these reference updates, 
the calculated equal-error threshold was augmented. This is perhaps 
not the correct strategy in light of the error-rate trend. 

From the third reference update, reference no. 4 and thereafter, the 
error rates decrease generally monotonically. It is reasonable to believe 
that a better choice of thresholds for the first three references could 
provide a strictly monotonically decreasing equal-error rate through 
the course of reference updates, say from 12 percent to 4 or 5 percent. 

To provide a reduced-error rate at initial stages, a larger set of 
initial samples, preferably recorded at different sessions, could be ob- 
tained. Another approach which may be more practical is to provide 
an additional decision category, that of “repeat’’ or ‘defer decision.” 
Thus, if the distance on any trial is within a specified fraction of the 
threshold, additional utterances are requested until a decision can be 
made. This strategy is useful because there is a considerable probability 
of acceptance on a trial after rejection on a previous trial. In this ex- 
periment, a tabulation of the frequency of acceptance following rejec- 
tion showed that 82 percent of the trials in which a customer was re- 
jected were followed by trials in which the customer was accepted. 
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(For well-adapted customer references, reference no. 5 and thereafter, 
this figure increases to 91 percent.) 

Another tabulation was carried out to assess the effect on error rate 
of withholding decision on trials for which |D/T — 1] < A, where 
D and T are the overall distance and threshold, respectively. For 
A = 0.05, the customer-reject rate fell from 8.8 percent to 7.4 percent, 
a 15-percent improvement, with decisions withheld on 4 percent of 
the trials. For the same value of A, a 22-percent improvement was 
noted for the accept-impostor rate. For A = 0.1, the customer-reject 
rate fell by 24 percent and the accept-impostor rate by 39 percent, 
with decisions withheld on approximately 7 percent of the trials. 

The system was shown to be capable of adapting to long-term varia- 
tion in speaking behavior by means of periodic updating of reference 
data and the use of an adaptive threshold. Regular and continuous 
usage of the system is probably necessary for such adaptation. Several 
customers had absences of two or more weeks during the course of the 
experiment but were not rejected with greater frequency on their 
return. However, the sample was too small to be conclusive. It seems 
likely, in fact, that many individuals will have difficulty being accepted 
with respect to old reference patterns after prolonged absences. 

The other principal ‘‘real-world”’ condition of interest was the effect 
of transmission and background noise over dialed-up telephone lines 
and the background noise at the calling location. Although this was 
not studied in any concentrated or organized manner, generally speak- 
ing, the system seemed quite tolerant of the transmission and back- 
ground conditions encountered. Over a nine-day period the following 
levels were monitored. The standard deviation of peak speech levels 
was approximately 8 dB. The background noise during recording was 
approximately 35 dB below the average peak speech level. Occasionally, 
background noise reached levels 10 dB higher than this average level. 
(A background noise level greater than —18 dB with reference to 
average peak speech level, or a peak speech level less than —12 dB 
with reference to average peak speech level, would prompt a request 
for a new recording.) Peak speech levels for female speakers were 
generally 2 or 3 dB below peak levels for males, while levels for calls 
originating from outside the local exchange were about 3 dB below 
the levels for local calls. (Approximately 25 percent of the total popu- 
lation called from outside the local exchange, generally via toll lines.) 
The only condition observed to be definitely detrimental was record- 
ing in the presence of pulse-like background noise, such as the kind 
originating from a typewriter or teletype console in close proximity 
to the telephone transmitter. One customer who habitually called 
under this condition is one of the outliers in Fig. 8. 
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Abnormal voice conditions were also monitored in an informal way. 
Mild upper respiratory infections were not observed to have any effect 
on the customer-rejection rate. However, a fairly severe case of 
laryngitis was observed for which the speaker could not provide an 
acceptable utterance because his voice “‘broke”’ at each attempt. The 
most severe voice condition problem observed was diplophonia. Di- 
plophonia is a condition associated with a husky or “raspy” voice 
quality. An inspection of the speech waveform for diplophonic in- 
dividuals reveals voiced speech intervals in which alternate pitch 
periods are more strongly correlated than adjacent pitch periods. 
Pitch analysis for such speakers is quite difficult. Two customers, one 
male and one female, were observed to be diplophonic. The female 
speaker provided a rejection rate of 44 percent. It seems clear that 
for any system making use of pitch analysis, diplophonic speakers 
should be identified and, if possible, pitch analysis modified or 
eliminated. 

It is useful to keep in mind the results of some previous supple- 
mentary experiments in the light of the present evaluation. One ex- 
periment assessed the ability of human listeners to perform a speaker- 
verification task.!° Using the same 8-customer, 32-casual-impostor 
speaker set used in the earliest evaluation of the automatic system, 
listeners performed at an average equal-error rate of 4 percent in a 
series of A-B comparison trials. That is, in 4 percent of the trials 
in which the speakers were the same, the listeners judged them to be 
different, and in about the same fraction of the trials in which the 
speakers were different, the listeners judged them the same. This per- 
formance level is somewhat better than what has been observed for 
this telephone evaluation, but not as good as previous laboratory 
evaluations. 

Another study! indicated that intensively trained professional 
mimics have considerably better chances for acceptance than casual 
impostors. With pitch-and-intensity analysis alone, mimic acceptance 
was found to be of the order of 15 percent when thresholds were set 
for a customer-reject rate of 1 percent. With the inclusion of pre- 
dictor coefficient analysis, mimic acceptance falls sharply to a tolerable 
4-percent level. Thus, pitch-and-intensity analysis by itself ina speaker- 
verification system may be rather vulnerable to the efforts of deter- 
mined mimics. 

The practical acceptability of a speaker-verification system depends 
on its expected error rate and its intended application. There are 
many noncritical screening applications in which the present error 
rate of 9 percent is acceptable. This error rate can be improved con- 
siderably (perhaps by a factor of 2 or more) with the use of individual 
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test phrases and multiple-phrase sequential strategies. For applica- 
tions with more stringent requirements, say error rates of the order 
of 1 percent and good protection against mimics, extended analysis 
techniques are necessary. As already mentioned, a laboratory evalua- 
tion which included predictor coefficient analysis satisfied these 
requirements. 


IV. SUMMARY 


The feasibility of operating the automatic speaker-verification sys- 
tem in the “‘real world’’ has been demonstrated. The system proved 
to be tolerant of many of the degraded and uncontrolled transmission 
and environment conditions which occur in the “real world’ when 
customers access the system from their own premises via dialed-up 
telephone lines. The error rate obtained for stable and established 
reference patterns is approximately 5 percent, which is quite accept- 
able considering the abbreviated analysis used in the experiment. The 
greatest weakness seems to lie in the establishment of adequate initial 
reference patterns. It is felt that at least a partial remedy for this 
difficulty can result from a collection of initial samples at more than 
one recording session and the use of a “deferred-decision”’ category if 
the distance is within a specified tolerance of the threshold. 
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Simulation of repeaters for an optical-fiberguide digital-communication 
system requires the calculation of the statistical error rate of the signal. 
The calculations of these error rates are difficult because of the non- 
gaussian nature of the noise in the optical detector. In this paper, statistical 
techniques are described which are useful for simulating tails of distribu- 
tions. In particular, the importance-sampling procedure is used to modify 
the probability densities of the input values in a way that makes simula- 
tion possible. Application of this procedure gives more accurate results 
in reasonable computer times. The method is applied to the calculation 
of the error rate of a fiberguide repeater. Realistic examples are simulated. 
Results compare favorably with experimental measurements. The number 
of samples needed for simulation was reduced by five to six orders of mag- 
nitude. An alternative numerical solution to the problem is also developed. 


I. INTRODUCTION 


Traditional computation of error rates of digital repeaters is often 
done by assuming that the noise is gaussian in nature. For a linear 
equalizer, the error rate at a given sampling time can then be easily 
computed if the signal and standard deviation of the noise are given.! 

The fiberguide repeater has essentially three kinds of noise sources. 
Two noise sources are located in the detector, while the third is a 
thermal gaussian noise generated by the input stage of the preamplifier 
(see Fig. 1). 

The detector, which is an avalanche photodiode is governed by two 
intrinsic noise processes: the primary photo current, which is a Poisson 
process, and the avalanche multiplication process which produces the 
gain in the diode. Neither of these processes is gaussian. 

The probability density functions (ppF) for the output current of 
the diode were derived by Personick?? and McIntyre.’ The expressions 
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Fig. 1—Optical repeater. 


describing these distributions are quite complex, and excessively time 
consuming for computer simulation of the repeater. An approximation 
to these probability densities was obtained by Webb, McIntyre, and 
Conradi (wmc).' The approximate ppr (referred to below as the 
WMC PDF) is, according to the authors, fairly accurate in the practical 
operating range, and is easily implemented on the computer. 

The wnc ppF for the current at the output of the diode is: 


1 —(xz — M)? 


p(x) = (+ EY eS) | 


(1) 








where 


M = (n.q/T)G = mean output current 
q = charge of an electron 
T = time interval 
G = average avalanche gain 
Ne = (nT /hQ)P.»t = number of primary electrons 
n = optical conversion efficiency 
Popt = average optical power in time T 
hQ = energy of a photon 
o? = n.G’F.(q/T)? = variance of the diode current 
A= VnF./ (F. 2m 1) 
F. = kerrG + [2 — (1/G)](1 — kere) = excess noise factor 
ketr = effective ionization ratio (assumed constant over range of 
G considered) 
kett and 7 are photodiode parameters. 


I 


Equation (1) is, more precisely, the ppF of the mean diode current 
over the time interval 7 for a given average optical input power, 
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Fig. 2—Probability density function for ONE and ZERO signals using wmc and 
gaussian approximation. 


P,pt, and a given average avalanche gain G. It corresponds to the sum 
over all primary electrons of the product of the primary electron ppF 
with the avalanche gain ppF. This equation, therefore, combines the 
effects of both detector noise sources. 

In Fig. 2 the solid lines represent current pp¥Fs for two signals: a 
ONE where the optical power is Pop, = —51 dBm and a ZERO 
where P.»t = —61 dBm. The other parameters are G = 80, y = 0.667, 
kere = 0.04, T = 22 ns. The broken lines represent error pprs of the 
same signals and same diode parameters where a normal approxima- 
tion to the distribution is assumed. The asymmetry of the wmc dis- 
tributions is clearly seen. 

The noise generated by the first stage of the amplifier has to be 
added to this detector-generated noise. 


Il. POSSIBLE SOLUTION TO THE PROBLEM 


To compute the error rate of the repeater, one has to know the 
signal probability distribution at the decision-making point, 1.e., at 
the output of the equalizer in Fig. 1. 
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The output of the equalizer with impulse response A(t) is 


y() = fea — adr, 


where (2) 
x(t) = X(t) + CTH (t) 
x,(t) = signal output of the avalanche diode 
Xru(t) = gaussian thermal noise generated by the amplifier. 


The problem is to find the distribution of y(t) for given distributions 
of x,(¢) and x7u(t) and impulse response h(t). 

Figure 3 shows the shapes of two distributions: Pi(y) is the cumula- 
tive distribution of y for a ONE signal transmitted, Po(y) is the prob- 
ability that the output current exceeds y for a ZERO signal trans- 
mitted. yo is the decision threshold of the comparator (Fig. 1). These 
distributions are 


y 
Py) = [ py/ONE)dy 
: (3) 
Poy) = |” p(y/ZERO)dy, 
y 
where p(y/ONE) and p(y/ZERO) are the respective conditional 
ppFs. It should be noted that, for practical laser sources, the ZERO 
input signal could be on the order of 10 percent of the ONE signal, 


LOG P(Y) 





Yort 


Fig. 3—Distributions for ZERO and ONE signals. 
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since the laser is not turned off completely. In addition, there is the 
dark current of the detector diode present. The error rate of a digital 
repeater is defined as 


P. = aPi(yo) + bPo(yo). (4) 


(The factors a and b are determined by the a priori distribution of 
the data. For unbiased data, a = b = 0.5.) 

One of the features that has to be determined is the optimum place- 
ment of the decision threshold yop, that is, where P, is minimum. 

A general analytical solution to this problem is not possible; there- 
fore, the problem has to be approximated numerically or simulated 
with the aid of a computer. Analysis of the system using a gaussian 
approximation to the detector noise sources was considered by S. D. 
Personick.® This analysis, however, predicts optimum decision thresh- 
olds considerably different from actual measured values.’? Another 
solution to this problem is to replace the exact solution by a simpler 
solvable formulation that gives an upper-bound solution. This was 
done by Personick*®:® using the Chernoff bounds. Finding upper 
bounds is a very effective technique since it gives a degree of confidence. 
Nevertheless, a method that gives more accurate results is needed 
to estimate how close the bounds are to the real solutions for this class 
of problems. 

The problem can be solved either by the numerical solution of eq. 
(8) or by statistical Monte Carlo simulation using importance sam- 
pling. Both solutions were attempted and are described below. Particu- 
lar emphasis was placed on the statistical simulation, since numerical 
solutions require excessive amounts of computer time. 


Ill. DISCRETE STATISTICAL ANALYSIS (DIRECT COMPUTATION OF 
OUTPUT PDF) 


The output of equalizer y(t) is given in eq. (2). Since we are interested 
in the output at the sampling time fo, it is not necessary to compute 
the convolution for every t. The value of to is selected to maximize 
y (to): 

yts) =f 2(r)h(o — 2)dr (5) 
since 
x(t) = x(t) + ron (t) 
y (to) i 2,(r)h(to — r)dr + / nou(r)h(te — r)dr 
Ys(to) + yru (to). (6). 
The problem is to find the probability distribution of y (to). 


I 
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Fig. 4—Input signal and impulse response of the equalizer. 


The probability distribution of yru(éo) is gaussian since try is a 
thermally generated noise. The standard deviation of yrx is 


4kT By 
OTH = Rex ? (7) 





where B,, is the effective noise bandwidth of h(t) and Rru is the 
thermal resistance. 
The effective noise bandwidth is defined as 


1 " . 
Be= rp [ |AUDas, (8) 


where H(f) is the one-sided Fourier transform of h(é) and |Ho| is 
the maximum absolute value of A(f). If ae impulse response is 
approximated by a raised cosine of duration 7, By = 3/4r. 

Typical waveforms for z,(t) and h(t) are shown in Fig. 4. The dis- 
crete approximation of the oa ys(éo) is: 


Ys = NV 7 Lnhny—n, (9) 


where 7’ is the duration of the signal and N is the number of samples 
used in the approximation. The samples zx, are independent and have 
a distribution of the form of eq. (1) for a time interval AT = T/N. 


3.1 Numerical computation of p(y) 


If the random variables x, are independent, then the density of 
their sum equals the convolution of their respective densities 


p( = sere | wo) wic)t 
Yel = ATR Thal (PM AT-hya) ? \ AT-hy_2 


= ( x.) |: (10) 
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Since the thermal noise is also additive and independent, the total 
density of the signal is 
p(y) = p(ys)*p(yru). (11) 


3.2 Fourier transform approach 


The simplest and fastest way to compute p(y) is to transform all 
partial densities using the FFT, multiply them all, and take an inverse 
FFT: 


PY) = TAP TT 7 Gor =e I || o( (ar ie) || 


-FLp(vra)]): (12) 


Unfortunately, the Fourier transform is computed as a sum of positive 
and negative numbers and, therefore, suffers from large numerical 
errors. Since we are interested in the tails of the distribution (error 
rate of 10-°), we conclude that the numerical errors will be excessive. 


3.3 Direct convolution 


Direct convolution produces less truncation error, because densities 
are positive functions. But direct convolution requires an excessive 
amount of arithmetic operations, e.g., if we have n distributions each 
having M points, the number of operations (multiplication and addi- 
tion) will be 

eo (13) 

Qn 

For M = 100 and n = 2, N = 104 which is a manageable number; 
but for M = 100 and n = 10, N = 3.5 X 10’ and the number of 
operations becomes excessive. A program using direct convolution 
was written for n = 2. A density for the average power of the input 
pulses was convolved with the density of the thermal noise. This can 
be regarded as a response of an integrate-and-dump equalizer (capaci- 
tive input). 

The results of this convolution are shown in Fig. 5. The solid lines 
are the distributions for a ONE and ZERO being transmitted. The 
average power of the ONE in this case is —51 dBm and of the ZERO 
—61 dBm, with T = 22 ns, Rey = 2.3 kilohms, and B, = 31 MHz, 
which corresponds to a ory = 14.68 nA. 

The decision threshold 6 is usually expressed in percent of the mean 
range of the output current as: 


Yo— Fo 


g= — 
¥i— Yo 


100 [percent], (14) 
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Fig. 5—-Distributions for ONE and ZERO signals with thermal noise; wmc and 
gaussian approximations. 


where 
Yi = mean output current for ONE 
Yo = mean output current for ZERO 
Yo = threshold current. 


As seen from the graph, the optimum decision threshold @ for this 
case is 42.3 percent and the minimum error rate is 9.86 & 107". 
These results were also observed in an experimental setup.’ The broken 
lines in Fig. 5 represent the same input but with the assumption that 
the probability density of the detector current is normal. As seen from 
the graph, the resultant threshold is considerably different from the 
one obtained with the wmc distributions. 


IV. STATISTICAL (MONTE CARLO) ANALYSIS 


Monte Carlo (mc) analysis has proven over the years to be a very 
effective method for analysis of systems whose parameters or inputs 
have known statistical distributions.° 
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However, straightforward application of Mc analysis for calculation 
of error rates is impractical. The reason is that the error rates of in- 
terest are in the order of 10—* to 10-°. To have some confidence in the 
results, 10 to 100 samples have to be in that region of low probability. 
This will require on the average 10! samples, which is impractical. 

To overcome this shortcoming, we introduced a modified sampling 
method called importance sampling.'° Conceptually, the idea is simple. 
If the regions which contribute to the result are known, modify the 
distributions in such a way that more samples are taken from the 
important regions. 

Let us consider the situation where the problem is to estimate the 
value of 


e= [- 9@)p@de, (15) 


where p(x) is a ppF. The importance-sampling procedure consists in 
introducing another por p*(z), which is preferable for sampling pur- 
poses. Then 





= foe) ES prez 


wr 
| 


i 


[ #@p* @ae. 
£ is then estimated from 

1X gt 

ans ug (zi), 
where 2; is picked from the ppF p*(z;). 


4.1 Defining the important regions for error-rate computation 


The problem is to define the regions where the errors come from. 
In Fig. 6a is shown the input waveform for a ONE, while Fig. 6b 
depicts the impulse response of the equalizer and 6c is the ppF of the 
diode current for one time interval AT. 

Since the input signal is always positive, it is obvious that for posi- 
tive h» most of the errors, that is, outputs that are below the threshold 
6, will come from region A in the ppr of Fig. 6c. Conversely, if hn is 
negative, most of the errors will come from region C. Exactly the 
reverse is true for transmitted ZERO signals, since the signal is con- 
sidered erroneous if it exceeds the threshold @. 

A very similar argument can be constructed for thermal noise, 
where the “important” regions are the positive tail for a ZERO trans- 
mitted and the negative tail for a ONE transmitted. 
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Fig. 6—(a) Optical input signal. (b) Impulse response of the equalizer. (c) PDF 
of average detector current for AT’. 


4.2 Biasing of the PDF 


To implement the importance-sampling technique, the distributions 
have to be modified so that the important regions are emphasized. 
The ppFs in the random-number generator that we use are stored in 
tabular form, as shown in Fig. 7a. The region in question can be 
biased by a constant multiplier as in Fig. 7b, where p* = pB for 
a <x Sb and where the bias B is constant. In Fig. 7c, the ppF is 
biased such that it is quasi-uniform in the region prescribed. The bias 
B; is constant in intervals between two adjacent points in the table. 
For each interval (x;_1, 2;), this bias is 


Bi (16) 
where 7; is the value of the ppF at x; and a is a coefficient specified by 
the user. The biasing and normalization is done automatically by the 
program. The output of the random-number generator consists of 
two values: a random number RV and its associated bias BV. BV is 
the bias of the interval that RV was picked from. 
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P(x} 





P(x) 





a Xi-1 Xj b 


Fig. 7—Probability density functions. (a) Unbiased. (b) Constant bias. (c) Bias 
for quasi-uniform PDF. 


4.3 Computing the histogram 

4.3.1 The usual histogram 

The standard output of a Monte Carlo analysis is the histogram. 
Usually the histogram is constructed by subdividing the output vari- 
able into M bins, Fig. 8. 

The statistical analysis program counts the number of outputs which 
fall into each bin. For a large number of samples (V — ~~), the histo- 
gram approximates the ppr. The probability p; that the output y will 
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Fig. 8—Histogram. 


be in the interval y;-1 S y S y;: is 


Di = Jim N? (17) 
where 
nm; = number of samples in bin 7 
N = total number of samples. 
The accuracy of estimation of p; can be established by" 
oi =p) 
= 2— = 2, /~__—_ 18 
lag Np a 


where +e is the 95.4-percent confidence limit of p;, p is the desired 
probability, and o; is the standard deviation of p;. 


4.3.2 Histogram for importance sampling 


Let us assume we have M independent random variables x, -:-, 
Lins a0, IM; 
y= f(x, ott) Bm) “s+, 0m), 


p* (2m) is the biased pp¥ of the random variable z,. How do we com- 
pute the unbiased histogram for y. 

If a sample xn; is selected from p*(%m), its probability to be selected 
is increased by the bias B,,;. Therefore, the weight of the sample zn; is 


— 
Bui ; 
Since y; is determined by M independent variables x, the probability 


of y:, to be selected from these regions, is increased by the product of 
all the biases of xm. Therefore, the weight of y;is w; = 1/([[ Bni). The 


Wmi = 


average weight of all the n; samples that fell into bin 7 is given by 
1 2 1 2 1 
Wavji = i Da Wi = — D ( ). (19) 
Jj 


44 n; =1\ I Bui 
m 
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The probability that the unbiased output will be y;-1 S y S y; is 
Nh; * 
Pj = Wavi Fy = WaviP;- (20) 


Since this method increases the probability p in eq. (18), the error 
e is reduced considerably. 


= }(1 =p") _ [Wav — Pp 


For error-rate computation p < 1, 


2 [wav 

eR = 2 : (22) 
VNp* Np 

and since pj is increased in the interesting areas by 1/wav;, the con- 
fidence limit is decreased on the average by 


y= 8 = Vows, (23) 


€C 








where 


ers = Monte Carlo analysis with importance sampling 


ec = conventional Monte Carlo analysis. 


4.4 Example 


A simple example to illustrate the principle of importance sampling 
is the problem of calculating the probability of obtaining a total of 
three when one tosses two dice. Each die has six faces labeled from 
one to six and each face has the probability p = 1/6 of being on top. 

The problem, of course, can be solved analytically. Any particular 
combination of the dice has the probability of 1/6-1/6 = 1/36 of 
occurring. Since there are two combinations which make three (one- 
two and two-one), the probability of getting a three in a random toss 
of the dice is 2/36 = 1/18. 

If the probabilities of a one and a two are biased so that they occur 
twice as often as usual, that is, p = 1/3 rather than 1/6, then the 
probability of getting a three will be four times as great as in the un- 
biased case. That is p = 2/9 instead of 1/18. The weight for this 
example is 

ws = 1/(B,B;3) = 1/4 


and the probability of obtaining a 3 is being approximated by 


Ai Syl 
Poa 


ale 
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The error in eq. (18) is reduced by slightly more than a factor of two, 
since for 


8.24 
=1/18 ¢=—= 
- / : VN 
and for 
3.74 
— 2 9, =—_lC 
e / . VN 


V. ERROR-RATE ANALYSIS 


The importance-sampling method described above was applied to 
the error-rate analysis of the fiberguide repeater with an avalanche 
diode detector. 

To demonstrate the analysis procedure, we will analyze one 
specific case. Let us assume that the input to the photodiode is as 
shown in Fig. 9. The ONE is a square pulse of —51 dBm power of 
22 ns duration. ZERO is 10 percent of the ONE power or —61 dBm. 
The avalanche diode data are 


average avalanche gain = 80.0 
ionization ratio k= 0.04 
conversion efficiency 7» = 0.667. 


The equalizer impulse response h(t) is as in Fig. 10 and the thermal 
input resistance is Rry = 2.3 kilohms. 

The pulses (ONE and ZERO) are subdivided into nine intervals 
[eq. (9) ], since this resolution is considered adequate for the impulse 
response in Fig. 10. Nine distributions* of the form in eq. (1) are com- 
puted for the ONE and ZERO pulses. 


Popt IN dBm 


-51 


t IN ns 
—| T=22ns ~— 


Fig. 9—Optical input. 
*In this case the distributions are identical since the pulse is rectangular. 
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h(t) IN OHM/s 


2.1014 


1.1014 


Fig. 10—Impulse response. 


The mean for ONE is M, = 400.38 nA 
The standard deviation for ONE is a, = 102.66 nA 
The mean for ZERO is M,= 40.03 nA 


The standard deviation for ZERO is oo = 32.46 nA 


The effective noise bandwidth [eq. (8)] for the thermal noise is 
B, = 31 MHz. The standard deviation is computed from eq. (7): 


OTH = 14.68 nA. 
From the above considerations, the output for N intervals is 


i> h aS h 24 
y= Hy At Non bop 2 nt TH, ( ) 


where the z, and xry are samples from the respective PDFs and h’s 
are discrete values obtained from h(t). The impulse response is normal- 
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ized to produce unity gain 





n=1 
n= > (25) 
De An 
n=1 
N * 
y* = ue hitn + Loe. 


In addition, the above case has been analyzed for N = 1. This cor- 
responds to the integrate-and-dump equalizer which was analyzed 
numerically in Section 3.3. 


5.1 Computation of distributions and biasing 


The calculations of distributions for the diode output are straight- 
forward using the wmc distributions of eq. (1). The biasing is a more 
complex procedure and sometimes requires a few trials. The idea is 
to get a reasonable amount of samples in the “important” region of 
the output. The ppr for ONE was computed for z, = 0 to 800 nA, 
and was biased according to eq. (16), 


B;= ar for 0 to 400 nA with a = 1. 

The ppr for ZERO was computed from 0 to 400 nA and was biased 
from 40 nA to 400 nA with an a = 0.1. The gaussian distribution for 
the thermal noise was tabulated for +7c. The biased regions were 
—7o to 0 for the ONE signal and 0 to 7o for the ZERO signal with 
an a = 0.5. The ppFs were stored in tabular form for 100 points. Since 
the ZERO signal has a ppF which requires more detail for small cur- 
rents (Fig. 2), the independent variable was spaced logarithmically. 
All others were linearly spaced. 

The distributions were computed by a discrete form of eq. (3): 


Cuml1(n) = > Hist (2) 


i=l 


(26) 
M . 
Cum0(n) = > Hist0(2), 
where Hist1(¢) and Hist0(z) are histograms of the outputs. 
To calculate the minimum error rate of eq. (4), Cuml(n) and 


Cum0(n) were smoothed and approximated by asymptotic func- 
tions. The minimum of the sum of these functions is the minimum error 
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rate, and the current where this minimum occurs is the optimum 
threshold. 


5.2 Smoothing of the cumulative distributions 

It was shown by 8. B. Weinstein” that for the class of ‘‘exponential- 
type” distributions the Pron prepayitly distribution has an asymptote 
in the tails of the form 


Pe) ~ exp[ -(2)]- (27) 


For »v = 2 the distribution is gaussian. 
To obtain better accuracy, this asymptotic function was expanded 
to 
P.(x) = exp{—exp[a(log x)? + blog x + c]} 
or (28) 
log (—log P.) = a(log x)? + blogz +c; 


i.e., the log log of P. is a second-order polynomial of log x. This equa- 
tion reduces to Weinstein’s form, eq. (27), for a = 0. 

The second-order polynomials were then least square fit to both 
Cum1 and Cum0. 


5.3 Minimum of error and optimum threshold 
The probability of error, eq. (4), is defined as 
PZ) = 2LP.() + P.(1)] 
= 3(exp{ —exp[ai(log I)? + bi (log I) + e1]}} 
+exp{—exp[ao(log I)? + bo(log I) + co]}), (29) 
where J is the decision threshold current. 
To find the minimum of P.(I), we compute Jop+ from 
dP..(1) 
dl I==T opt 
PAL rnin = Peroni): 


Equation (30) is solved using the Newton-Raphson method. The 
threshold @ is usually expressed in percent as in eq. (14). 


= 0. (30) 


VI. RESULTS 

Figure 11 shows a histogram for the ZERO signal for 10,000 samples. 
The abscissa is scaled logarithmically. In the tables on the right side 
are indicated the ordinate values for the bins, the number of samples 
in each bin, and the distribution for the bin. Below the histogram are the 
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NO. OF CUM. 0 
BIN (nA) SAMPLES DISTRIB. 


1 5.000E +01 85 2.232E—03 

2 6.000E+01 111 6.067E—04 

3 7.000E+01 165 4.128E—04 

4 8.000E+01 168 1.667E—04 

5 9.000E+01 216 9.279E—05 

6 1.000E+02 277 2.061E—05 

7 1.100E+02 293 1.573E—05 

8 1.200E+02 331 9.328E—05 

9 1.300E+02 339 6.318E—06 
10 1.400E+02 383 5.809E—05 
1 1.500E+02 427 7.079E—08 
12 1.600E+02 390 4.162E—08 
13 1.700E+02 456 1.349E—08 
14 1.800E+02 454 7.077E—09 
15 1.900E+02 458 1.984E—09 
16 2.000E+02 455 3.315E—10 
17 2.100E+02 440 1.114E—10 
18 2.200E+02 433 2.208E—11 
19 2.300E+02 444 5.407E—12 
20 2.400E+02 426 3.579E—-13 
21 2.500E+02 375 1.461E-13 
22 2.600E+02 380 6.224E-14 
23 2.700E+02 319 1.003E—14 
24 2.800E+02 307 8.280E-15 
25 2.900E+02 252 9.154E-16 
26 3.000E+02 208 7.867E—16 
27 3.100E+02 178 9.662E—17 
28 3.200E+02 194 1.057E-17 
29 3.300E+02 156 1.557E—19 
30 3,400E+02 142 1.030E—19 
31 3.500E+02 113 8.583E—-20 
32 3.600E+02 98 7.897E—21 
33 3.700E+02 76 4.744E—22 
34 3.800E+02 49 2.739E—23 
35 3,900E+02 54 2.839E—24 
36 4.000E+02 42 1.640E—24 
37 4.100E+02 29 8.339E—25 
38 4.200E+02 18 8.271E—25 
39 4.300E+02 20 5,535E—27 
40 4.400E+02 13 2.754E—29 


—25 —20 -15 —10 —5 
DENSITY (LOG SCALE) 


COEF. OF LEAST SQ. FIT Ag = 0.9954 A, = —1.2324 A> = 0.53019 
NUMBER OF SAMPLES = 10,000 


Fig. 11—Histogram for ZERO signal. 


coefficients of the fitted polynomial RW) = ao + ai log J + az(log I)?, 
where R(Z) = logio [—logio Peo (J) ]. 

Figure 12 shows the smoothed distributions for the above example 
for one section (N = 1) and 1000 samples. The stars represent nu- 
merical solutions for the same input parameters, as described in Section 
III. As seen from the data, if biasing is chosen rightly, a comparatively 
small sample size is adequate to simulate the extreme tails of the ppF. 

Figure 13 shows the smoothed distributions when h(t) is represented 
by nine sections (V = 9). The broken lines represent the distributions 
for N = 1. The sample size was 10,000 in both cases. As seen from 
the graph in this example, an integrate-and-dump equalizer analysis 
can be used as a first-order approximation to an equalizer design with 
an impulse response similar to the h(t) shown in Fig. 10. 
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The minimum error-rate and optimum threshold were 


N=9 Pemin = 3.03 X 107% 6 = 43.61 percent 
N=1 Penn = 6.6 X10-" 6 = 42.68 percent. 


Figure 14 again shows the distribution of the above example (N = 9). 
The circled points are results obtained using Chernoff bound approxi- 
mations as obtained by J. H. Bobsin® using Personick’s? method. As 
seen from the graph, the optimal threshold predicted by the Chernoff 
approximation is identical with the one obtained by the statistical 
analysis. The minimum error rate shown by the Chernoff bound is 
about two orders of magnitude larger from the statistically simulated 
value. However, the error rate obtained by the Chernoff approximation 
is an upper bound. Therefore, if the result falls within the acceptable 
design specifications, no further analysis is necessary. If, on the other 
hand, the predicted error rate is marginal, a more accurate analysis 
technique such as the one described in this paper is advisable. 

Figure 15 depicts curves for the minimum error rate Pemin and 
optimum threshold as a function of the standard deviation of the 


ag 


—5 


IMPORTANCE 
SAMPLING 


LOG P (I) 


O NUMERICAL 
CALCULATION 





0 50 100 150 200 250 300 
| IN NANOAMPERES 


Fig. 12—Distributions for ONE and ZERO signals. 
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LOG P (1) 


MONTE CARLO 

(1 SECTION) 
— MONTE CARLO 

(9 SECTIONS) 





0 50 100 150 200 250 300 350 400 
| IN NANOAMPERES 


0 25 35 40 45 50 55 60 75 100 
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Fig. 13—Distributions for ONE and ZERO signals using Monte Carlo method for 
1 and 9 sections. 
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(9 SECTIONS) 


LOG P (I) 


—--— CHERNOFF BOUNDS 





0 50 100 150 200 250 300 350 400 
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25 35 40 45 50 55 60 75 
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Fig. 14—Distributions for ONE and ZERO signals using Monte Carlo method and 
Chernoff bounds. 


764 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1976 


E1 — ERROR RATE FOR 1 SECTION 
E9 — ERROR RATE FOR 9 SECTIONS 
61 — THRESHOLD FOR 1 SECTION 

69 — THRESHOLD FOR 9 SECTIONS 


LOG Pemin 


Oop IN PERCENT 





0 5 10 15 20 25 30 
O74, IN NANOAMPERES 


Fig. 15—Minimum error rate and optimum threshold for 1 and 9 sections as a 
function of thermal noise. 


thermal noise ory. For this example, 


23.32 
= nA ], 
omm = [nA] 


where Rey is the input impedance of the equalizer in kilohms. 

The solid lines represent Pe min and Oop for N = 9. The broken lines 
depict the same case analyzed with N = 1. 

As seen from the graphs, the optimum threshold for both types of 
receivers 1s between 42 and 46 percent. Interestingly, the threshold has 
a minimum at o = 25 nA (Rey = 0.87 kilohm). 

The integrate-and-dump analysis has a lower error rate then the 
nine-section equalizer analysis in the normal operating range of the 
repeater. 

Throughout this analysis, the avalanche gain was assumed to be 
constant. In a more realistic design, the avalanche gain can be adjusted 
to optimize the error rate. This will result in a higher avalanche gain 
for higher thermal noise. 

In this example, the input signals and the equalizer impulse re- 
sponse were chosen so that there is no intersymbol interference. This 
was done to simplify the example. However, the analysis procedure 
is not limited to this case. 
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Vil. CONCLUSIONS 


The statistical importance-sampling technique was used to calculate 
error rates of optical repeaters. The application of this technique sub- 
stantially reduced the number of samples needed for simulation of the 
repeater and thus made computer simulation practical. 

The importance-sampling technique is very effective and has general 
applicability in areas where tails of distributions have to be computed 
or particular sensitive areas of the statistical distribution have to be 
emphasized. 

This method cannot, however, be used automatically ; a good prior 
knowledge of the behavior of the system is necessary. To get an 
accurate statistical simulation of a section of an output, all the im- 
portant regions of input ppFs which may contribute to this section 
have to be emphasized. 
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A Phased, Optical, Coupler-Pair Switch 


By V. RAMASWAMY and R. D. STANDLEY 


(Manuscript received December 15, 1975) 


We analyze an optical switch consisting of two directional couplers 
with an intermediate phase-compensating network. The conditional re- 
quirements on coupler length, phase mismatch between guides in each 
coupler, and phase compensation between couplers are established such 
that two switching states exist. It 1s shown that, to switch between states, 


(t) The couplers must be 3 dB or stronger. 
(22) Kach coupler must be ‘“‘adjusted’” to become a hybrid by control- 
ling the phase mismatch between guides. 
(212) The differential phase @ must be initially set to turn on one state. 
(tv) An additional increment of x in the phase shift causes switching to 
the second state. 


It is estimated that the controlling voltage required to activate the bipolar 
switch is of the order of 1 to 2 V for a 3-um guide with 38-mm-long elec- 
trodes in Ti-diffused LiNbO3. We are currently in the process of fabricat- 
ang such devices. 


1. INTRODUCTION 


A number of methods!~‘ have been suggested for realizing optical 
switches in integrated optical format; in many cases, the schemes are 
extensions of the microwave art.>® Recently, Papuchon et al.” and 
Campbell et al.® obtained switching by perturbing the phase syn- 
chronism in indiffused strip-guide directional couplers of about one 
coupling period in length. Kogelnik et al.® have proposed a switch 
consisting of tandem couplers with adjustable phase mismatch be- 
tween guides in each coupler. This approach results in relaxed dimen- 
sional tolerances over previously described coupler switches. However, 
the adjustment of phase mismatch in both couplers is required to 
switch between states. 

In this paper, we propose and analyze an optical switch consisting 
of directional coupler pairs with an intermediate phase-compensating 
network. We show that the switching between states can be achieved 
by a m change in the intermediate phase shift providing that each 
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coupler has been initially adjusted for hybrid performance, i.e., for 
3-dB coupling. Prior to adjustment, the coupling strengths are required 
to be 3 dB or stronger, implying coupler lengths of 4 to 14 basic 
coupling periods. 


Il. PRELIMINARY ANALYSIS 


We assume directional couplers with uniform coupling having trans- 
fer matrices of the form” 


[ei]= [a 4-][e). ® 


where |A|? + |B|? = 1 (see Fig. 1). For coupler 1, the matrix coeffi- 
cients relating the field amplitudes H;, are given by" 





Ai = cosa, — ; 1. sin a, 2 
1 1 J Vy +1 1 ( ) 
. sma, 
B= jj, 3) 
oni Vyi t 1 ( 
and, for coupler 2, 
As = COS az — ee 4 
2 a2 a ar | Oe, ( ) 
‘ sin ae 
Be = j ——= (5) 
EET! 
where 
a12 = V¥i2 + 1 wha (6) 
2Lo 
V1.2 = ABs 2b (7) 
T 


AGB.,2 = propagation constant difference between the guides in 
couplers 1 and 2, respectively, 


the interaction length, 


lie 


Lo 


transfer length for total power transfer for y = 0. 


If we connect the two couplers in series via a differential phase-shifting 
network (Fig. 1), then, for unit excitation of port 1 of the first coupler, 
we obtain for the output signals 


Es = A,Ae aes BiB.e*, (8) 
E, = —A,Bt — AtBie%, (9) 
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Fig. 1—Two couplers with adjustable Ag, connected in series with an intermediate 
phase shifter. 


where ¢ is the differential phase. The two states of the switch are 


State I: E;=0, |#.| =1 
State II: |£;| = 1, E, = 0. 


Thus, for State I, the condition is 


| 


AiA2 





OL = 
e/ BiB, ) (10) 
i.e., 
AAo| _ 
BB,| = 1 (11) 








This simply implies that, if the couplers have coupling coefficients 
that are the inverse of each other in magnitude, then a unique ¢ = ¢1 
exists that permits the realization of State I. 

For State II, the condition is 





ABs 
$2 => — — 
ere? A3B;’ (12) 
1.e., 
AB3| _ 
AtB*| ~ 1. (13) 








This implies that, if the magnitudes of the coupling coefficients are 
identical, then there exists a ¢ = ¢: that permits the realization of 
State IT. 

It follows, then, that with the parameters of the one coupler speci- 
fied, there exist values of ¢ that permit the realization of States 
I or II for each allowed value of yz and a2. To realize both states by 
adjustment of ¢ alone, hybrid couplers are necessary, i.e., |Ai] = 
[Ae] = [Bil = |Bo| = 1/v2. 

Low-loss strip waveguides made by diffusion of metals in electro- 
optic materials’-"“ are available to facilitate the fabrication of these 
switches. However, the difficulty in forming coupler pairs is the 
maintenance of uniform guide width, separation, and refractive in- 


PHASED OPTICAL SWITCH 769 


dex differential.* Thus, realization of hybrid pairs is rather difficult. 
The above analysis suggests that, given couplers with a power-coupling 
coefficient greater than 3 dB, States I or II can be achieved by appro- 
priately altering the phase mismatch in the couplers (electro-optically 
or otherwise) and introducing a specific differential phase shift ¢ be- 
tween couplers. 

Thus, our requirements are 


(¢) Each coupler must initially have coupling coefficient |B1\, 


|B.| = 1/v2. 

(22) Direct-current biasing can then be provided in each coupler to 
alter the phase mismatch in order to achieve |B,| = |B! 
=1/v2. 


(7iz) A differential phase shift ¢1 is then introduced between 
couplers for State I to cxist, and then additional phase shift of 
T, Viz., 62 = 1 + (2m + 1)z, causes the switch to go to State 
II. Thus, activation of either state can be achieved by the 
addition or removal of the control voltage that causes a phase 
shift of « between couplers. 


Section III provides a detailed analysis and establishes the initial 
conditions to realize the switch. 


Ill. DETAILED ANALYSIS 


In this section, we derive the equations governing the states of the 
switch in terms of the coupler parameters and the differential phase 
shift. 

Substituting eqs. (2) through (5) in (10), we obtain 


eid = E + cos? a1 3 + cos? =a ei {tan —%Y1 tan = | 


sin? a4 sin? a (1 + yi)! 
| 72 tan os || 
+ tan G+ \} (14) 
The phase shift ¢; required to null port 5 is 
—yi tana — 72 tan ae 
= tan- Yi 1 t —1 | — | ; 15 
nies | atopy | Tl ae! ae 
and the additional condition established in eq. (14) is 
vyitcostea, sitar. 
sinta, 2 + cos? ae (16) 


Equation (16) relates the coupler parameters such that we have 


* The effect of dimensional tolerances on individual coupler performance has been 
treated by several authors (Refs. 15-18). The effect of such tolerances on interference- 
switch performance has been considered by Shelton (Ref. 19). 
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couplers whose ratio of the power-coupling coefficients are the inverse 
of each other in magnitude, i.e., |A1|/|Bi| = |Bz|/|Az2|. Inspection 
of (15) and (16) reveals that, for a given value of y1, a1, there are 
many combinations of y2, az, and ¢; that will meet the inverse con- 
dition requirement and therefore realize State I. 

Equations (2) through (5) can be substituted in eq. (12) to yield 
the following conditions for the phase shift ¢2 required to null port 6: 


es 23 | =r tan ay _, | —Ye2 tan ae 
$2 = 7 + tan | ee ve + tan | eens ans (17) 


Using eq. (15), (17) can be written as 
gg = dita, (18) 
and the additional condition is 


vit cos? a1 _ 3 + cos* a (19) 
sin? ay sin? a» 


Equation (19) illustrates the need for the ratio of the coupling coefh- 
cients to be identical, i.c., |Ai|/|Bil = |A2|/|B2|. Again, we have 
many combinations yz, a2, and ¢2 that will meet this requirement to 
realize State II for a given value of y1, a;. Although (18) is written in 
terms of ¢1, note that the values of y1, a; and the associated values 
of y2, a2, and @¢2 are distinctly different from those needed to realize 
State I. 

Obviously, if we set |B:| = |Be| = 1/v2, i.e., if the couplers be- 
come hybrids, then eqs. (16) and (19) are simultaneously satisfied and 
both states can be simultaneously realized by simply changing ¢ by 
adding zw. In other words, the de biasing in each coupler is so chosen 
that (y1, a1) and (ye, a2) satisfy 


| sin Q1,2| ot 1 . 
(ie + DIV pe 
For illustrative purposes, we consider only the fundamental solution 
of eq. (20) for positive y, which is plotted in Fig. 2. The plot shows the 
normalized length of the coupler as a function of the normalized 
phase mismatch y to achieve hybrid performance. In practice, the 
couplers as fabricated are not usually exact hybrids. For any phase 
mismatch between guides, the coupler lengths have to be within } 
to 1} transfer lengths (Zo). In practice, this can be easily achieved. 
To illustrate the adjustment for hybrid performance, let us consider 
two couplers with the following initial specifications: 


Coupler x: yz = 0.24, I/Lo = 1.2, 
Coupler y: yy = 1.1, I/Lo = 1.1. 


PHASED OPTICAL SWITCH 771 


| 
CIN 






\ 


FORy>1 
HYBRID PERFORMANCE CANNOT BE 
ACHIEVED 


0.6F 


PHASE MISMATCH vy 


Fig. 2—Required coupler length (//Zo) to obtain hybrid operation as a function of 
the phase mismatch parameter y of the coupler. 


As seen from Fig. 2, biasing the two couplers to change y’s by Ay. 
and — Ay,, respectively, results in hybrid performance, with y., = 0.58, 
vy = 0.7. Then, initializing the phase shifter to a value ¢1 given by 
eq. (15), the switch can then be operated by changing the value of 
¢1 by adding z. 

To obtain the phase change gui, let us further assume the two couplers 
are identical in performance, i.e., y1 = y2, a1 = a2. (This assumption 
is easily satisfied, especially with fabrication processes using electron- 
beam lithography with analog scan or step-and-repeat techniques.) 
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Fig. 3—¢1 and ¢2 required to switch on State I and II, respectively, assuming 
identical coupler pairs and both being adjusted to provide hybrid operation. 


For this case, then, 


gi: = 2 tan7 | waa | (21) 


and 
g2 — G1 = T. (22) 
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Equations (21) and (22) are plotted in Fig. 3. Assuming Ti-diffused 
LiNbO; strip waveguides,” the control voltage required to produce 
a phase shift of 7 is estimated to be 1 to 2 V for a 3-um guide width 
and L = 3mm. The special case with ¢ = 0 is treated in the appendix. 


IV. CONCLUSIONS 


Equations have been established governing the states of a switching 
network consisting of two dephased, directional couplers interconnected 
via a differential phase network. The results imply that, in practice, 
one would strive to fabricate coupler pairs where power-coupling 
coefficients are greater in magnitude than 3 dB and, by appropriate 
bias, adjust each for 3-dB coupling. Two unique values of differential 
phase will then exist such that both switching states can be realized. 
The voltage required to switch between states is estimated to be of the 
order of 1 or 2 V. 


APPENDIX 


If we set ¢: = 0, then we have the equations governing the switch 
described in detail by Kogelnik et al.?! It can be shown that, in addi- 
tion to requirements of eq. (16) (State I), we must have 


2 sin? a; 
sa Bie, a eh 23 
1 7 + cos? a1 (28) 

and 
af-y2 4 
vi + cost ay 

tan = 2 —.—_ - 24 
= sin? ay (24) 


These simply imply that, given y; and a, y2, and a, are uniquely fixed. 
This is because ¢ = 0 is a unique state and the second coupler must 
satisfy both magnitude and phase conditions. 

The graphical data presented by Kogelnik are quite useful in that 
the magnitude requirement of eqs. (16) and (19) is implied in those 
results.”! 
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Mode Mixing With Reduced Losses in 
Parabolic-Index Fibers 


By D. MARCUSE 
(Manuscript received January 13, 1976) 


We present design criteria for the construction of a modified parabolic- 
index fiber with intentional mode coupling. Mode coupling serves the 
purpose of reducing multimode pulse dispersion and ts accomplished by 
introducing carefully designed index fluctuations into the fiber core or 
by controlled ‘‘random’”’ bends of the fiber axis. Radiation losses due to 
mode coupling can be minimized by terminating the parabolic-index 
fiber core in an abrupt index discontinuity. The additional modes intro- 
duced by this step must be filtered out by periodic mode filters that consist 
of parabolic-index fiber sections without the refractive-index step. 


Il. INTRODUCTION 


In multimode optical fibers, light pulses are carried by many modes. 
Such multimode operation limits the information-carrying capacity of 
an optical fiber communications system because of pulse spreading. 
At the fiber input, all modes receive part of the energy of the light 
pulse simultaneously. However, at the fiber output, the light pulse is 
stretched out in time because each mode (or mode group) travels at 
a different group velocity. The length of the stretched pulse is pro- 
portional to the fiber length. The amount of pulse spreading in a multi- 
mode fiber depends on its construction. The dependence of the group 
velocity of the modes on the mode labels is influenced by the distribu- 
tion of the refractive-index profile of the fiber in radial direction.! 
Step-index fibers* typically exhibit more pulse spreading than fibers 
with graded refractive-index distributions. Fibers whose index profiles 
show a parabolic (or square law) dependence on the radial coordinate 
have the property that all modes have nearly the same group velocity 
so that pulse spreading on parabolic-index fibers is nearly minimized.! 

Our remarks so far apply to multimode transmission in perfect fibers 
where each mode propagates independently of all the other modes. 


*Step-index fibers have a core with constant refractive index and a lower index 
cladding. 
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Inhomogeneities of the axial distribution of the refractive index or 
geometrical imperfections of the fiber geometry tend to couple the 
modes among each other. Mode coupling has the undesirable con- 
sequence that some of the power is coupled to unguided radiation 
modes resulting in radiation losses. On the other hand, coupling of 
the guided modes results in a continuous interchange of power be- 
tween fast and slow modes so that a light pulse that is distributed 
over all the modes assumes a new shape. Instead of spreading according 
to the different group velocities of the independent modes, the light 
power carried by coupled modes is forced to travel at an average 
velocity and exhibits a narrower width, which spreads only propor- 
tionally to the square root of the fiber length.? Thus, mode coupling 
can be intentionally introduced to improve the pulse performance of 
multimode fibers.4 However, this technique for improving the pulse 
performance must be used with great care to avoid an unacceptable 
increase in the power loss of the fiber. 

To understand how mode coupling can be tailored to minimize losses, 
it is necessary to consider the coupling process in more detail. Each 
mode has a characteristic propagation constant By. The label M is 
used to identify the mode. A unique mode designation requires that 
M consists actually of two symbols representing the radial and azi- 
muthal mode numbers. For simplicity, we combine the double label 
in the single symbol M. As mentioned before, mode coupling is pro- 
vided by some deviation of the fiber from its perfect geometry and 
composition.® We use a function f(z) to describe the axial dependence 
of the deviation of the refractive-index distribution or of the core- 
radius deviation from their nominal, perfect values. In addition to 
this function of the length coordinate z, we need its Fourier transform 
which we define as 


F(6) = lim lz i : fle)eide|- (1) 


Coupling between two modes labeled M and N is mediated by a 
particular Fourier component of f(z) according to the law? 


Bu — Bn = 8. (2) 


It is thus clear that two modes remain uncoupled if F(@) = 0 for the 
particular 6 value required by (2). Furthermore, if the differences be- 
tween the propagation constants of neighboring guided modes depend 
on the mode number, F'(6) # 0 is required over a certain range of 
6 values if all guided modes are to be coupled. Using these rules and 
certain “selection rules,” it has been demonstrated‘ that most guided 
modes of a step-index fiber can be coupled with very little radiation 
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loss if the Fourier spectrum (1) is limited to a carefully selected range, 
such that Fourier components exist for coupling the guided modes 
according to (2), but that coupling between guided and radiation 
modes [also obeying the law (2) ] is prevented. 

Mode coupling in parabolic-index fibers requires only a narrow spec- 
trum of spatial frequencies @ because the differences between propaga- 
tion constants of neighboring modes are almost independent of the 
mode label. This feature of the modes of the parabolic-index fiber 
causes a problem since it makes it harder to discriminate between 
coupling among guided modes and coupling from guided to radiation 
modes. With an appropriate selection rule, coupling among guided 
modes of the step-index fiber has the property that the differences (2) 
increase with increasing mode number.‘ Cutting the Fourier spectrum 
(1) off at a maximum spatial frequency 0 = @max thus stops mode 
coupling at a given mode number, so that coupling between guided 
modes of lower order and modes with the highest mode numbers is 
prevented. Since only the highest-order modes are near (in mode 
number space) radiation modes, coupling between guided and radia- 
tion modes is thus avoided. Because of the almost constant differences 
between propagation constants of neighboring modes, this strategy 
fails in the parabolic-index fiber. 

We show in this paper that we may couple the guided modes of a 
parabolic-index fiber and still avoid radiation losses by using a modified 
parabolic-index profile. An ideal parabolic-refractive-index profile has 


the form 
n=n[t—(ZYal Ir] < @. (3) 


r is the radial coordinate and A/a? determines the gradient of the 
index profile. The true parabolic-index profile cannot be realized be- 
cause the refractive index of ordinary solid materials cannot be less 
than unity. Typical parabolic-index fibers have index profiles of the 


form 
2 
m[1- (5) 4] 
_ a 
~ |e = no(1 — A) a 
1 b 


lA 


a 
; (4) 


oO, 


n 


IIA IIA 
/\ IA 


The cladding region a S |r| Sb is usually so thick that, mathe- 
matically, we may assume b > o. The guided modes do not carry 
significant amounts of power inside of the cladding region so that 
they behave almost as if they were guided by the ideal index distribu- 
tion (3). Modes with significant amount of power in the cladding region 
are no longer guided but belong to the continuous spectrum of radia- 
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tion modes. By providing a narrow Fourier spectrum of spatial fre- 

quencies for the purpose of coupling the guided modes, we necessarily 

couple the highest-order guided modes to radiation and lose power. 
The situation is changed if we modify the index profile to the follow- 


ing form: 
Yr 2 
_ {no [1 ~ (zYa| || 
n= a 


Ne, < no(1 — A) a < |r| 


IIA 


(5) 


ioe) 


IA 


(b = © was assumed for simplicity). The index profiles of (8), (4), 
and (5) are shown in Figs. la, lb, and 1c. We divide the guided modes 
of the index profile (5), shown in Fig. lc, into two classes. There are 
modes whose field distributions are essentially limited to the region 
0 < |r| Sa. These modes have negligibly small field intensities in 
the region |r| = a and behave as though they belonged to the idealized 
medium defined by (3). Because modes of this kind are essentially 
the modes of the parabolic-index medium, we call them P-modes. 
There are, in addition, modes of order higher than the P-modes whose 
field distributions reach strongly into the region near r = a. These 
modes are guided by the index discontinuity at r = a and behave 
similar to the modes in a step-index fiber. For this reason we call them 
S-modes. The differences (2) between neighboring P-modes are nearly 
identical, while the differences between S-modes are much larger and, 
if an additional selection rule is introduced, increase with increasing 
mode number. If a narrow band of spatial Fourier components is pro- 
vided to couple the P-modes among each other, S-modes will remain 
uncoupled from P-modes and also remain uncoupled among each other. 
Thus, we have achieved coupling among the P-modes and have im- 
proved their pulse performance. However, if S-modes were allowed to 
reach the detector, the pulse performance of the fiber would be de- 
graded very seriously because of the different group velocities of 
P-modes and S-modes and the large group velocity spread among the 
S-modes. It is thus necessary to suppress S-modes before they reach 
the detector. This can be done easily by adding at the end of our fiber 
with index distribution (5) (or Fig. 1c) another fiber section that does 
not allow S-modes to propagate. A fiber with an index profile accord- 
ing to (3), or Fig. 1b, has this property. 

The strategy just outlined would work if there were truly no coupling 
between P- and S-modes. Some residual coupling is, however, unavoid- 
able because of imperfect fiber tolerances. A small amount of power 
will always be coupled from P-modes to S-modes causing a ‘‘noise”’ 
background to reach the detector. This unwanted noise can be reduced 
by installing mode filters periodically along the fiber. As mentioned 
above, mode filters for S-modes consist simply of fiber sections with an 
index profile according to (4), or Fig. 1b. If we construct a fiber whose 
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Fig. 1—(a) Ideal parabolic-index profile of eq. (3). (b) Truncated parabolic-index 
profile of eq. (4). This profile is used for mode filters. (¢) Modified parabolic-index 
profile of eq. (5). This profile is used for the fiber with intentional mode coupling. 


index profile is given by (5), or Fig. 1c, for most of its length, but 
which is changed to assume the form of (4) shown in Fig. 1b for rela- 
tively short sections periodically interspersed with the rest of the fiber, 
we obtain a fiber guide with mode filters for S-modes. Mode filters must 
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not be spaced too closely in order to avoid excessive losses. The addi- 
tional loss occurs because we cannot avoid coupling between P-modes 
and a small group of S-modes that lie immediately adjacent (in mode 
number space) to P-modes. Along the boundary between P- and §- 
modes, mode spacing cannot be controlled so that we must assume 
that the intentionally introduced strong coupling mechanism will 
couple P-modes to their immediate S-mode neighbors along their 
common mode boundary. The mode filters strip away all S-modes, 
thus causing a small amount of loss of power that has been coupled 
to the S-mode group near the mode boundary. Design criteria for 
optimum mode filter spacings will be given in this paper. 

If intentional mode coupling is achieved by introducing index 
fluctuations into the fiber core, the mode filters can be incorporated 
into the fiber by the same manufacturing process that was used to 
produce the intentional fiber “‘imperfections.’”’ On the mode filter 
sections, no strong coupling will be provided to avoid additional losses. 
We shall show later that coupling can also be provided by small bends 
of the fiber axis. 

In the following chapters, we provide the necessary information to 
explain the mechanism and give design criteria for a modified parabolic- 
index fiber with mode coupling. 


Il. MODE SPECTRUM OF THE MODIFIED PARABOLIC-INDEX FIBER 


We have defined P-modes and S-modes in the introduction. P-modes 
have field distributions that have (negligibly) small values at the core 
boundary r = a. Their properties are nearly identical with the modes 
of the ideal, infinitely extended square-law medium. We may associate 
rays with these modes. The rays corresponding to P-modes spiral 
around the fiber axis in helical paths. Axial rays cross the fiber axis 
and move out to a turning point where the tangent to their path is 
parallel to the fiber axis. Spiraling rays encounter two turning points, 
one near the fiber axis and the other one at larger radii.!:* An approxi- 
mate field description utilizes the fact that modes in weakly guiding 
fibers are nearly linearly polarized’ (at least suitable superpositions 
of exact fiber modes are nearly linearly polarized). Their dominant 
transverse electric field component can be approximately expressed by 
the WKB approximation in the form! ® 


peak cos (v = eee (6) 
({L(k)r? — B?}r? — v?)# 
with 


ip es i ({[n(r)ke — g2}22 — »®)Par. (7) 
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A is an amplitude coefficient, the integer v is the azimuthal mode 
number, @ the propagation constant, and 


= w (€ojt0)?, (8) 


the propagation constant of plane waves in free space. The inner 
turning point r: is the smaller root of the equation 


{[n(@r)&P — Br? — = 0. (9) 


Since the denominator of the expression in (6) vanishes, the WKB 
approximation seems to fail at the turning points. A more careful 
analysis (not used here) is capable of bridging the gap and extending 
the validity of the wxB solution across the turning points. The require- 
ment of continuity of the field solutions across the turning points 
leads to the condition® ® 


[75 Ula @e? — Be = ar = p+ 0) 


The upper turning point r, is the larger root of eq. (9), the integer p 
is the radial-mode number. Equation (10) defines the propagation 
constant 6 and is called the eigenvalue equation. For P-modes we 
have re < a and obtain by substitution of (5) or (3) into (10) the 
expression 


3 
B= {nae - noe (2A)#(2p +v+1)7 - (11) 


S-modes reach the core boundary. It is well known that an abrupt 
index change at the core boundary forces the electric field intensity 
to assume very low values at r = a. For all modes, with the exception 
of modes very close to cutoff, it is permissable to approximate the 
actual value of the electric field by H = 0 at r = a. This condition in 
conjunction with (6) and (7) leads to the following eigenvalue equation 
for S-modes, 


[5 UeOeP = yt — r= p45 CY) 


Substitution of (5) and integration leads to 


[ (koa)? — 2(noka)?A — v2} 
—y | arctan (er Stee) & 4 


(koa)? : (2noka)?A — (koa)? 7 
2noka (2d)? | aresin (ex = SES OaT ) my 5 | 


=(2p + 3)m (13) 
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with 
Ke = nek? — B?. (14) 


The guided and radiation modes of the modified parabolic-index 
fiber can be displayed in mode-number space. Each mode is charac- 
terized by two integers, the azimuthal mode number » and the radial 
mode number p. Mode-number space displays the values of v and p 
in the plane shown in Fig. 2. For P-modes we can introduce a com- 
pound mode number 

M=2p+v+1 (15) 


and express the propagation constant (11) in the form, 
3 
B= | (nok)? — 2™% (aay (16) 


Modes with constant values of M have the same propagation con- 
stants and le on diagonal dotted lines, two of which are shown in 
Fig. 2. Cutoff for the S-modes is defined by the condition 


B = nok (17) 
or 
Ko = (mg — n3)}k. (18) 


When we substitute (18) into (13), we obtain an implicit equation 
for the mode boundary between S-modes and radiation modes shown 
in Fig. 2 as the solid line. The boundary between P-modes and S-modes 
is defined by the condition in (16) written in the form 


(koa)? = 2noka(2A)?*M (19) 
and from 
(koa)? = 2(noka)?A + v?. (20) 


Equation (20) expresses the requirement that the transition from 
P-modes to S-modes occurs when the turning point falls on the core 
boundary and is obtained by combining (5), (9), (14), and the con- 
dition r = r. = a. By eliminating «x ,a from (19) and (20), we obtain 
a functional relation between v and p which was used to calculate the 
dash dotted mode boundary shown in Fig. 2. This figure was computed 
with the help of the following parameters, 


Nh, = 1.5 
ka = 150 et) 
A = 0.0098 


These values lead to the following value for the refractive index at the 
core boundary: mo(1 — A) = 1.515. 
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Fig. 2—Mode number space. The solid line is the boundary between S-modes and 
radiation modes. The dash-dotted line is the boundary between S-modes and P-modes. 
The dotted lines labeled M = 9 and M = 18 are curves of constant compound mode 
number M. M = 18 is used as approximate mode boundary. 


For the purpose of designing a coupling mechanism for the P-modes, 
we need to know the differences between the propagation constants 
of the modes. If we let » change by é6y and p by 6p, we compute from 
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the eigenvalue equation for the P-modes in the form (19), 


Oko oe 6M 
mo 2M oe 
with 
5M = 26p + 5S». (23) 
For the S-modes, we find from (13) 
a 2 42 
a(dp) + 4(dp) | arctan (er —*) coe 
oko _ vs 2 (24) 
ae 1 2 a2 
(2p + $)r + vj arctan so) 45 —S 
vs 2 
with the abbreviation 
S = [ (koa)? — 2(noka)?A — v?]}. (25) 
The difference of the propagation constants follows from (14): 
2 
58 = — = au (26) 
Ko 


We have plotted values for | 68a| in Fig. 3. These values were computed 
from (22) and (26) for P-modes and from (13), (24), and (26) for 
S-modes using the numbers in (21). We also assumed that 6M = +1. 
For P-modes, 6M = +1 clearly yields the least separation between 
guided modes. (6M = 0 would lead to 68 = 0 and is excluded.) Just 
as in the case of the step-index fiber, we introduce a selection rule by 
properly designing the coupling mechanism (see Section III below). 
The selection rule is 


év = +1. (27) 
We may now achieve a transition between neighboring modes with 
|6M| = 1 in two different ways. We can either use 
éy = +1, bp = 0 (28) 
or 
5yv= +1, dp= ¥1, (29) 


where the upper or lower set of signs belong together respectively. All 
other combinations lead to bigger values of |6M@| and larger differ- 
ences for the separation between the S-modes. 

Figure 3 shows that the values of | 68a| for P-modes are very nearly 
independent of the mode number M. By using a Fourier spectrum (1) 
with spatial frequencies in the range 6a = 0.14 to 0.1415, we couple 
all neighboring P-modes with |6M| = +1. The hatched area in Fig. 3 
labeled ‘‘S-modes”’ indicates the range of | 68a| values that occurs for 
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Fig. 3—The differences of propagation constants of adjacent modes with |sM| = 1 
as functions of the compound mode number M. 


a given value of M. It is apparent how widely the differences between 
S-modes with |6M| = +1 vary. Typically, the combination (28) leads 
to smaller 68 differences. However, even the smallest spacing between 
adjacent S-modes is so much larger than the corresponding spacing 
between P-modes that it should be relatively easy to design the 
coupling mechanism so that S-modes are not intentionally coupled 
among each other. The differences between adjacent P-modes and 
S-modes along the dash-dotted mode boundary in Fig. 2 cannot be 
calculated from our simplified theory, so that we must assume that 
P-modes, lying along the mode boundary, may be coupled to their 
S-mode neighbors. 

We have now proved that P-modes can be coupled among each other 
by a spatial Fourier spectrum of very narrow width and that it is 
relatively easy to discriminate between P-mode and S-mode coupling. 
It remains to discuss mode coupling, study the degradation of pulse 
performance that results from coupling a few S-modes along the mode 
boundary to P-modes, and to discuss the performance of the mode 
filters. 
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lll. MODE COUPLING AND PULSE SPREADING 


Mode coupling is provided by implanting a refractive-index per- 
turbation into the fiber core.? Instead of the index distribution n given 
in (5), we now use a perturbed refractive-index distribution n so that 
we have 
n? — n? = < f(z) cos ¢ lr| Sa. (30) 
The dependence of (80) on cos ¢ imposes the selection rule (27). 
Without such a selection rule we could not uncouple P-modes and 
S-modes. The linear r dependence indicated in (30) is quite arbitrary. 
It is necessary that (30) vanishes at r = 0 in order to have a well de- 
fined value of the function at this point. The linear r dependence is 
not only the simplest function that vanishes at r = 0 but also one for 
which the coupling coefficient can be evaluated. Other r dependent 
functions could be used as factors in (30). The function f(z) is assumed 
to be a random function with the narrow Fourier spectrum (1) that 
assures coupling among the P-modes but prevents coupling between 
P-modes and S-modes. 

The system of coupled P-modes can be described by coupled-power 
equations.2 The coupling coefficients are obtained by solving overlap 
integrals including the function (80) multiplied with the electric field 
functions of the two modes whose coupling is to be evaluated. The 
P-modes are described very well by the Laguerre-Gaussian field solu- 
tions of the ideal parabolic-index distribution (8).9"° We obtain for 
the power-coupling coefficients? : 


yopw—1.p = K(p + ») (31) 
op y—1,p41 = K(p + 1) (32) 
op w¥1,p = K(p eae 1) (33) 
hyp +1 ,p-1 =, Kp (34) 
with 
k 
K (35) 


= T6n8a(2A)! (F?(6)). 


The symbol ( ) indicates that an ensemble average has been taken. 
The argument of the Fourier spectrum function (1) is the appropriate 
difference of the propagation.constants of the coupled modes according 
to (2). [See eq. (75). ] 

It is our intention to solve the coupled-power equations for the 
P-modes with the coupling coefficients given by (31) through (35). 
The coupling mechanism postulated by the refractive-index perturba- 
tion (80) couples a given mode (», p) to its neighbors (vy + 1, p), 
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(vy ~1,p+1), (» — 1, p), and (v + 1, p — 1). In the first two cases, 
the M number is increased by one while it is decreased by one in the 
third and fourth case. The coupled-power equations thus assume the 
form? 


OP, 1 oP, 
ae Pp + a ar + aPy» => expel aa _ P55) 


hin sip pct pha P55) + | ne © ene ard Py) 
Pip vio P iat _ Pia) (36) 


We have indicated by our notation that the group velocity vj, depends 
only on the compound mode number JM defined by (15) and not on 
the individual values of » and p. In addition, all modes are assumed 
to suffer the same loss a caused by absorption and random scattering 
processes. 

To be able to solve equation system (86), we assume that modes 
with equal values of the compound mode number M carry equal 
amounts of power. It is not obvious that this must be the case, but 
we may argue that modes with equal mode number M have identical 
propagation constants and, thus, are coupled among each other by 
the zero-spatial-frequency component of F(@). The zero-frequency 
components of random distortions tend to be very large. Consider, 
for example, the departures from perfect straightness of the fiber axis. 
No fiber in actual use is ever perfectly straight. In fact, its slow varia- 
tions tend to be particularly large so that its power spectrum peaks 
at zero spatial frequencies. This fact has been observed whenever 
power spectra of fiber distortions have been measured.”:!3 Conse- 
quently, we assume that modes with equal values of M are strongly 
coupled and, hence, equally populated by a coupling process that is 
not explicitly taken into account. It is only incorporated into the 
analysis by assuming that modes with the same values of M carry 
equal amounts of power. Making this assumption and substituting 
(31) through (35) into (36) allows us to write, 











FP, = Py 
Pap = Peta = Pati (37) 
Pity = Pyeapsi = Pye 
and 
oP 1 oP 
a — —* + ePy = KU(M + 1) (Panui — Pu) 
2 Un ot 


—M(Py — Pu-s)]. (88) 


If the number of guided modes is very large, we may regard them as a 
quasi-continuum and treat M approximately as a continuous vari- 
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able.'* This assumption allows us to write the equation system (38) 
as a partial differential equation,’ 








OPy 1 0Py_ _ a 74, aPu 
a2 ‘om ob oP + Kop | MSH | (39) 


Using (16), we can approximate the inverse group velocity in the 
following way 


1 _ 1968 _ n , 
ue eek ee (40) 
with the abbreviation 
A 
~ noc (ka)? a) 


c is the velocity of light in vacuum. We have indicated in Section I 
that a group of S-modes along the mode boundary between P- and 
S-modes are coupled to the P-modes. These modes have different 
group velocities. We incorporate this mode group into our theory by 
assuming that the modes with the highest value of M, M = N, have 
a group velocity that differs from the law given by (40) for P-modes. 
Thus, we use 


+ 4 WrM? + pN23(M — N)]- (42) 
Um Cc 


The parameter p is the relative amount by which the inverse velocity 
difference 1/v,, — no/c of the mode group with M = N differs from 
the normal value. If the S-modes followed the law for the P-mode 
group velocity, we would have for the mode group with M = N 
<= 04 WH. (43) 
UN Cc 
Equation (42) states instead that the highest-mode group has inverse 
group velocity 
S = 4 (1 + p) WN. (44) 
UN Cc 
We solve the time-dependent coupled-power equations in the usual 
way.” First, the time-independent problem is solved. The only differ- 
ence from the normal procedure is that, in our present case, we do 
not assume that the highest-mode group is depleted by a loss process. 
Thus, instead of requiring Py = 0, we assume that power outflow 
stops at M = N and require correspondingly 


OP x 


aM =0 at M=N. (45) 
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The time-independent problem, thus, has the solution 
Pap = en *# Sg, LOU MN) 
i=0 VN Jo(us) 


The parameters w; are the roots of the Bessel function Ji(z) = 0 with 
uo defined as zero; thus, 


—aie, (46) 


ui = 0, 3.832, 7.016, 10.123, ---. (47) 
The eigenvalues are 
2 
o=K (48) 


with K of (35). The coefficients c; determine the initial power distribu- 
tion at z = 0. For large values of z, P reaches its steady-state solution 


Co — LZ 
P M VN é , (49) 
which indicates that all modes carry equal amounts of power. 

The time-dependent, steady-state solution is known to represent a 
pulse with gaussian shape.?* Its full width (in time) between 1/e 
points is given by second-order perturbation theory? and can be calcu- 
lated to assume the following form: 


r= 1)(1+%): (50) 


For p = 0, we obtain the pulse width of P-modes that are completely 
independent of the S-modes, 





= Qn d/2A sian sie 
Po = VL Tea TUR) TE 


Lis the length of the fiber that has been traversed by the pulse. The 
pulse width spreads only proportional to the square root of the dis- 
tance traveled. Equation (50) shows that the group of S-modes that 
is coupled to the P-modes widens the pulse more if its group velocity 
departs more from that of the P-modes—that is, with increasing values 
of p. Its influence is reduced with increasing values of the maximum 
mode number N. 

A comment needs to be made regarding our analysis. We have 
calculated the impulse-response width (50) by treating the mode 
boundary between P- and S-modes as being parallel to the lines 
M = const. Figure 2 shows that this is not strictly true. In addition, it 
is not true that all S-modes along the mode boundary have equal 
group velocities. Our result must thus be regarded as an estimate and 


(51) 
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we must use an average group velocity for the S-modes. For M = N, 
we use the value that is obtained by setting v = N/1.5 in (19) and 
(20) and obtain, 

N = [0.573nokav2A Jint- (52) 


The subscript “‘int’’ attached to the bracket indicates that the nearest 
integer to the number in brackets must be taken. The approximate 
mode boundary (52) is shown in Fig. 2 as the dotted line labeled 
M = 18. 


IV. THE INFLUENCE OF RESIDUAL COUPLING BETWEEN P- AND S-MODES 


If the intentionally introduced coupling were the only mechanism 
by which the guided modes interact, we would have no further problem. 
However, residual coupling between all the modes is unavoidable so 
that we must consider the problem of power coupling from P-modes 
to S-modes and vice versa. 

If the P-modes carry light pulses, some of their power will couple to 
S-modes. Since S-modes travel with different group velocities, the 
power they receive from P-modes spreads out and forms an almost 
continuous background of noise power. Some of this power is coupled 
back into P-modes so that this noise background reaches the detector 
even if we filter out the S-modes before they reach the fiber end. 
Periodically spaced mode filters will reduce this noise problem, but 
filters increase the overall losses of the system so that it is necessary 
to reach a compromise between excess loss due to mode filters and 
undesirable noise caused by unintentional S-mode coupling. 

We treat the noise problem in two stages. First, we consider the 
power-coupling process from P-modes to S-modes, and in the next 
step we allow this power to couple back to the P-modes. 

Consider the equation system? 


OQ n 1 dQ, 
0z ae ot 








=), fal & (53) 


These equations describe only coupling from P-modes to S-modes. 
Coupling in the reverse direction is ignored as are losses. For weak 
coupling in short, low-loss fiber sections, this equation system is a 
reasonable approximation to the complete coupled power equation. 
H represents the residual, undesired coupling mechanism; the sum- 
mation extends over all P-modes and the indexing system is simplified 
by using one symbol for the complete set of mode labels. By using a 
double Fourier integral transformation, it is possible to derive the 
following solution of (58), 


Q.= 5 (2 A nH) , =e Jew (' Te) —~ af (a) |, (54) 
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which is based on the assumptions that the pulses carried by the P- 
modes travel with a uniform velocity v (these pulses are strongly 
coupled by the intentional coupling mechanism) and that their shape 
is gaussian,” 


P, = A, exp |- poe | (55) 


with energy content 


| oan ae ge (56) 


and that the S-modes are slower than the P-modes, 
Wn <»v forall n. (57) 


The shape of the pulse carried by one S-mode is shown in Fig. 4 
as the solid line. The shape of the P-mode pulses is indicated by the 
dotted line; the relative height of the two pulses is of no significance. 
The width of the S-pulse is 


6z = t(v — wn) = (V — Wn) —. (58) 


Because the P-pulses are spaced as closely as possible to maximize 
the information rate, we may safely assume that the S-pulses overlap 
after a distance z that is short compared to the total length of the 
fiber. The number of S-pulses that overlap at any given point is, on 
average, 

_ 8 _ v7 Wn (59) 


POWER 


Spe ESE _- P—PULSE 





z 


Fig. 4—The solid line is the shape of the S-pulse that results from coupling of 
P-mode power whose shape is indicated by the dotted line. The relative height of 
the two pulses is of no significance. 
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The average power carried by a given S-mode is, thus, 
Q., = MQ, = =0> pe A yyLiy. (60) 


[ The difference of the two error functions in the bracket in (54) assumes 
the value 2 over most of the region where its value does not vanish. 
Equation (60) can also be obtained directly from (53). ] 

We have now determined the average power in a given S-mode due 
to unintentional power coupling from the P-modes. We have also 
convinced ourselves that the S-modes do not carry pulses, because of 
extensive pulse overlap, but carry a quasi-continuous background noise 
signal. 

The amount of noise power that is coupled back from S-modes into 
P-modes is obtained from the time-independent coupled-power 
equations? 

= 





This equation is solved by ones integration. Summing over all modes 
and substituting (60) yields the following expression for the total 
noise power that is carried by all the P-modes, 


oy 2: 
No= XP, Ho RCO VT Mn mE). (62) 


Greek summation symbols are used to indicate summation over 
P-modes while Latin symbols express summation over S-modes. 

We assume that the fiber length between mode filters is 2 = d. At 
the end of the fiber of length ZL, we collect the noise contribution from 
L/d fiber sections between mode filters. The total noise is, thus, 


v= 34 ELL HH lly (63) 


So far we have ignored absorption and scattering losses. If all modes 
suffer identical losses, we only need to multiply (63) by the total loss 
that a signal suffers in traveling through the fiber. However, since our 
objective is to derive an expression for the signal-to-noise ratio, the 
absorption loss drops out in the end because signal and noise suffer 
identical losses. The average signal power at the end of the fiber is 


1 
S=F x Hy. (64) 
The signal-to-noise ratio at the fiber end is, therefore, 
S&S 2D eae (65) 
Ni vTLd PB Se : A 2 OF 
vu n 
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If the unintentional coupling coefficients Him were known, we could 
calculate the signal-to-noise ratio from (65). Equation (65) tells us 
that the signal-to-noise ratio improves as the spacing D between pulses 
is increased or the pulse width T is decreased. The signal-to-noise 
ratio deteriorates with increasing unintentional coupling strength, in 
fact this decrease is inversely proportional to the square of the coupling 
strength. The signal-to-noise ratio decreases linearly with total fiber 
length Z and with increasing fiber length d between mode filters. It 
is thus advantageous to space the mode filters as closely as the addi- 
tional loss that filters introduce will allow. 

For practical purposes, it will be necessary to obtain information 
about the unintentional mode coupling by making signal-to-noise 
measurements on a representative fiber sample. The signal-to-noise 
ratio of a fiber sample of length Z, follows from (65) if we set L = d 
= [,. We may use pulses of length 7’, for this measurement (7, or T 
are measured at the end of each fiber) and use a pulse spacing D, for 
the calibration measurement. This allows us to express the signal-to- 
noise ratio of a fiber with mode filters in terms of the measured signal- 
to-noise ratio (s/n), of a representative fiber sample, 


S _ DT.L; 
Ni ~ D.TLd TLd (s/N)c. (66) 


V. LOSSES CAUSED BY MODE FILTERS 


We explained in the introduction that mode filters are required to 
suppress the buildup of power in the S-modes. A mode filter is simply 
a section of fiber—without the mechanism for intentional mode 
coupling for P-modes—whose refractive index profile is modified from 
the shape of Fig. le (used for most of the fiber) to the shape of Fig. Lb. 

Mode filters introduce additional losses because the strong coupling 
mechanism provided for the P-modes couples these modes to a group 
of S-modes immediately adjacent to the common mode boundary. 
Let us first consider how many modes there are for a given compound 
mode number. If we count the number of combinations of vy and p that 
lead to a fixed value of M defined by (15), we find the following 
expression, 


Nw =| 5" | +1. (67) 





The subscript ‘“‘int’’ indicates in this case that the integer smaller than 
the number in brackets must be taken. For simplicity we use the 
approximation 


New: (68) 
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As a further approximation, we replace the mode boundary by the 
straight line labeled M = 18 in Fig. 2. The total number of P-modes 
is now approximately 

N M 

aes mS 4N(N + 1). (69) 

Because of the strong intentional coupling, each P-mode and each 

mode in a group of Nw S-modes along the common mode boundary 

carries the same amount of energy. The energy reaching the cth fiber 

section is /'7,. The total energy loss on the cth fiber section just after 

the mode filter is, thus, 

Ny 

Epis — Ere = — | = 

ats in Noy + Nw 

The symbol a indicates the losses caused by absorption and random 
scattering. We solve the difference equation and use the relation 


eted 4 (1 — ct | Er,. (70) 





od =2 (71) 
to obtain with the help of (68) (with M = N) and (69), 
a aa 1 —2az 
an ee (Gres) | (72) 


At the end of the fiber, we have a total of L/d mode filter sections. 
The excess loss caused by the mode filters is, thus, 


Err 2aL — N + ves) — p~2asL 
Ee.” eS = gntash, (73) 


This allows us to define the approximate power-loss coefficient per unit 
length for the excess filter-loss penalty, 


Bee BO 
Nd nokav2A d 
Equation (52) was used to obtain the right-hand side of this expres- 


sion. Short fiber sections of length d between mode filters (in other 
words, more mode filters) thus increase the excess loss. 


Zar == (74) 


VI. DISCUSSION AND NUMERICAL EXAMPLES 


In the introduction, we have outlined our strategy for reducing pulse 
spreading in parabolic-index fibers. We are now in a position to offer 
specific design criteria for our approach. The most important aspect 
of the fiber design consists in incorporating index fluctuations into the 
fiber core that obey the relation (30). As pointed out earlier, the r 
dependence of the index perturbations is not particularly important. 
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The linear dependence was chosen for convenience. The function f(z) 
is a nearly periodic, random function that must contain spatial fre- 
quencies in the range 


_ 2a 
a 
with 69 = 1.15 6A 


[This formula follows from (2) and (26) with the help of (16), (18), 
(22), with 6M = 1, and (52) ]. The parameter A is defined by (8). 
For the numerical example given by (21), we find that 6a (a is the 
fiber-core radius) covers the range from 0.14 to 0.1415. To gain insight 
into the meaning of these numbers, we assume that the vacuum wave- 
length of the light transmitted through our fiber is \ = 1 ym. This 
means that ka = 150 leads to a core radius of 23.87 um. The spatial 
period length of the almost periodic random function must thus be 
A = 1.07 mm. 

Next, we must decide what amount of index fluctuation is required 
to achieve a desired reduction in the width of the uncoupled pulse 
(pulse carried by uncoupled modes). In the absence of coupling, the 
pulse width is given by 


V24 2A 
a aa (75) 


_ nol as 

Puc => “oC A’. (76) 
The pulse width in the presence of mode coupling must be smaller 
than this value, otherwise steady state has not been reached and (50) 
and (51) are not applicable. First, let us consider how much harm 
results from coupling between P-modes and the unwanted group of 
S-modes along the common mode boundary. Using the values in (21), 
we find from (52) that N = 18. Let us (arbitrarily) assume that p = 1; 
we then find from (50) that the pulse width of the coupled modes is 
28 percent larger than the width of coupled P-modes in the absence 
of S-mode coupling. Ignoring the slight pulse-broadening effect of 
residual S-mode coupling, we define an improvement factor R as 


To 


ee Fe 


(77) 





It is desirable to make R as small as possible and only values with 
R <1 are meaningful. From (51), (52), and (76) we get 


14inai 0.83 
(2@) (cr) 


The numerical value on the right-hand side was obtained with the 


R= (78) 
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help of (21). To achieve, for example, R = 0.1 would require 


| Zo] = 3.3. 


To understand what this result means in terms of the amplitude of 
the refractive-index fluctuations, let us consider that the function f(z) 
appearing in (80) has the form 


f(@) = Asin [2 + y(2)], (79) 


where y(z) is a random phase function with correlation length Dorr. 
The power spectrum of this function defined by (1) is* 


sin? | @ — Q) _ 


(F?(0)) = A? (0 — Q)2Deor : 


(80) 


Instead of the correlation length Deo, we may introduce the width 
50 of the spectral band by the relation 


4a 





od 7 Deorr ey 
With Q = 6, we thus obtain from (80) and (81) 
L } L\} 
| z (F(6)) | =A (4,3) = oS: (82) 


(The number on the right-hand side pertains to our example.) This 
important relation shows how the amplitude A of the refractive-index 
fluctuations is related to the spectral bandwidth 66, the core radius, 
and the length of fiber over which a certain improvement factor R 
(in our case R = 0.1) is to be achieved. For our example, #a ranges 
from 0.14 to 0.1415 so that we have aéé = 0.0015. Assuming L = 1 
km, and using a = 23.87 um, we have from (82) 


A=1.1 X 10-5. (83) 


The refractive-index fluctuations follow from (30) and (79) if we 
assume that n — n <1 [1 is the perfect index distribution (5) ]: 


ee ay ey ee 
n-n= oa a sin [6z + W(z)_] cos ¢. (84) 


It is apparent from (83) that very slight index fluctuations are very 
effective for intentional mode coupling. Much more substantial pulse- 
length shortening than R = 0.1 should, thus, be easily achievable. 
The random phase fluctuations with correlation length Deorr can be 
produced by keeping the phase y(z) of (79) constant over a fixed 
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distance and introducing a random phase jump periodically in length 
intervals Deorr. The relation between D.or and the desired bandwidth 
of the random function is given by (81) and (75). 

Not much has been said about the problem of actually implementing 
the refractive-index fluctuations prescribed by (30) or (84). In prin- 
ciple, it is possible to introduce the refractive-index fluctuation during 
the process of preform fabrication, since it 1s necessary to take special 
care to produce the parabolic-index profile in any case. It may thus be 
possible to produce the intentional deviations from the perfect para- 
bolic-index profile in the core by programming the process of chemical 
vapor deposition, or whatever process is used to control the refractive- 
index profile. However, there is a much simpler way of realizing a 
refractive-index fluctuation of the kind required by (80). Let us write 
the index profile inside the fiber core in cartesian coordinates: 


n= no{1—- SEE a], (85) 


a 





If we displace the index profile in x direction from its symmetric posi- 
tion, we may make the substitution 


LIL — YG. 


Assuming the displacement g to be small, we obtain instead of (85) 


n=n-+ 2no— A (86) 


with n once more given by (85). By transforming the cartesian co- 
ordinates to cylindrical coordinates using 

x= rcos 9, (87) 
we obtain 


n? — n? XY 2no(n — n) = - (ania 3) cos ¢. (88) 


Comparison of (30) and (88) allows us to make the following associa- 
tion: 


I(@) = Anig 2. (89) 


The actual displacement of the fiber axis from perfect straightness can 
be expressed as 


g = Bsin [Qz + Y(z)], (90) 
with B being the displacement amplitude. Comparison of (79), (89), 
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and (90) allows us to use the relation 


aA 
B= Tea un 


For our example values (21), we thus obtain from (83) 
B= 2.9 X 1073 um. (92) 


This discussion shows that we may introduce the desired perturbation 
by actually bending the fiber axis in the form indicated by (90) with 
the very small amplitude B. Such bending could be accomplished by 
surrounding the fiber with a suitably strained plastic jacket. 

The abrupt refractive-index discontinuity of magnitude no(1 — A) 
— 2 is necessary to insulate the P-modes from the radiation modes by 
creating a buffer region of S-modes. The amount of this index step is 
actually quite arbitrary. However, it is clear from the general formula 
(65) for the signal-to-noise ratio that the sum over n in the denominator 
of this expression is larger when it consists of more terms, that is, if 
there are more S-modes. To keep the signal-to-noise ratio as small as 
possible requires keeping the number of S-modes small, which is 
achieved by keeping the abrupt index discontinuity small. What is the 
minimum index step that is needed for effective isolation of the 
P-modes? We can estimate a minimum index step by the following 
consideration. The lowest possible value of the propagation constant 
of P-modes is given by 

Bm = no(1 — A)dk. (93) 


The largest propagation constant of the radiation modes (actually, 
what we call radiation modes become cladding modes in a fiber with 
finite cladding thickness) is given by 


Br _ Nok. (94) 


Adjacent groups of P-modes with the same value M of the compound 
mode number are spaced a distance @ (in 8-space) apart. A sufficiently 
wide buffer zone of S-modes is required to isolate the P-modes from 
radiation modes. Thus, we require that the P-mode group with propa- 
gation constant (93) is separated (in 6-space) by 36 from the radiation 
modes. This requirement leads us to the desired condition for the 
minimum height of the index step 


7] 2 
no(l — A) — my = SP = aL (95) 


For the values listed in (21), we find mo(1 — A) — ne = 0.008. 
It remains to consider the design of the mode filters. Mode filters 
consist of fiber sections with the refractive index distribution of Fig. 
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Fig. 5—Refractive-index profile of a fiber providing continuous mode filtering for 
S-modes by electromagnetic tunneling. 


1b and serve the purpose of stripping off S-modes. The length of the 
filter sections depends on the S-mode losses of the filters. Since S-modes 
are leaky waves on the filter sections, their losses are very high, so 
that filters 10 cm to 1 m long should be sufficient. We have seen that 
the signal-to-noise ratio (66) increases as the distance d between mode 
filters decreases. The only limitation on the length of d or the number 
of mode filters L/d is the additional loss caused by the filters. If we 
decide to tolerate an additional filter loss of 1 dB/km (2a; = 0.23 km“) 
for the entire fiber, we find from (74) with the numbers of (21) 


d = 0.47 km. (96) 


If the filter spacing is reduced below this distance, the filter-loss con- 
tributions increase above 1 dB/km. 

The signal-to-noise ratio that is caused by unintentional coupling 
between P-modes and S-modes can be computed from (66) if a calibra- 
tion measurement has been made. In principle, it would be possible 
to compute this value from (65); however, such a calculation would 
require a detailed knowledge of the unintentional power-coupling 
coefficients H,,. In the absence of accurate information of unintentional 
fiber imperfections, such a calculation is not possible. 

Instead of using discrete mode filters spaced at certain intervals it 
is also possible to design the fiber so that mode filtering is achieved 
continuously over its entire length. The index profile of such a fiber 
is shown in Fig. 5. 
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We show formally how the application of boundary layer techniques to 
a well-stirred, concentrated, multicomponent electrolyte yields essentially 
the same results as those for a dilute, binary electrolyte. In particular, 
the coupled integral equations for the reactant concentration and the 
electric potential, derived by J. L. Blue for a dilute, binary electrolyte, 
remain valid in the general case with modified values of certain combina- 
tions of physical constants appearing in the equations. 


1. INTRODUCTION 


In this paper, we use the method of matched asymptotic expansions 
to derive boundary layer equations for a well-stirred, concentrated, 
multicomponent electrolyte containing a single ion reacting at the 
electrodes. We have two main objectives. First, we wish to show how 
the boundary layer approximation, up to now applied mainly to dilute, 
binary electrolytes or to the case of excess supporting electrolyte,! 
may be systematically generalized. Second, and more important, we 
show how the generalization leads to equations almost identical with 
those valid for a dilute, binary electrolyte, for which J. L. Blue? has 
derived an elegant and efficient method of solution. 

In the general case, the singular perturbation in the boundary layer 
leads to a set of coupled convective diffusion (cp) equations for the 
concentrations or, in other words, a vector cb-equation in place of the 
single, scalar cp-equation for a dilute, binary electrolyte. The vector 
cp-equation is diagonalized by introducing the eigenvalues and eigen- 
vectors of a reduced diffusivity matrix. The resulting uncoupled, scalar 
cb-equations can then be solved by Blue’s method. 

The perturbation yields linear cp-equations in the electrolyte, but 
the strongly nonlinear dependence of the electrode current density on 
the electrolyte potential and the reactant ion concentration is retained. 
In fact, the first approximation is simply the so-called ‘secondary 
current distribution,” the potential satisfying Laplace’s equation in 
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the electrolyte, but satisfying nonlinear boundary conditions at the 
electrodes. Furthermore, if, as is often the case, the current-carrying 
capacity of the electrode surface reaction is large compared with that 
due to ohmic conduction, the potential is approximately constant on 
the electrodes, giving the so-called ‘‘primary current distribution.”’ 

We begin in Sections II and III with a semiqualitative discussion 
of some of the approximations involved. Sections IV through IX 
contain the detailed analysis. The reader not interested in the details 
may skip to Sections X and XI, which contain very brief summaries 
of the formal perturbation solution and the equivalent boundary 
integral solution developed by J. L. Blue. 


ll. THE FIRST APPROXIMATION 


A first approximation to the cathode current density serves as a 
convenient introduction to the general concentrated, n-component 
electrolyte. In a well-stirred electroplating bath, the ion concentrations 
are almost constant, except in thin diffusion boundary layers surround- 
ing the electrodes. If the electrode current densities are bounded, so 
are the concentration gradients in the boundary layers. Thus, the 
changes in concentrations across the boundary layers are small when 
the boundary layer thicknesses are small enough. To a first approxi- 
mation then, the concentrations are constant throughout the elec- 
trolyte, current being carried entirely by ohmic conduction. In the 
electrolyte, the dimensionless potential y (equal to F®/RT, for Fara- 
day constant F, gas constant R, and absolute temperature 7’) satisfies 
Laplace’s equation 

Veo = 0. 
On the electrodes, the current densities are determined by surface 
reactions producing cations (e.g., Cut*) at the anode and depositing 
them at the cathode, the reaction rates being functions of the local 
concentrations and potentials. The potential y then satisfies the 
boundary condition 


— Ja; on the anode, 
k©do/dn = kOn-Ve = Aes on the cathode, 
0, otherwise, 


where n is the outward normal, k® is the dimensionless bulk con- 
ductivity, and, for e = a, ¢, 


je = je Lexp (ae y — gel) — exp(—Be| ¢ — vel) ]. 


The dimensionless ion exchange current density j;, here constant, 
but in general a function of the concentrations, and the exponents a, 
8. characterize the surface reaction rates. The quantity ¢ — ¢- is the 
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so-called ‘‘electrode surface overpotential’” for electrolyte potential 
and electrode potential gz. 

When j- is large, | ¢ — ¢-| must be small or, to a first approxi- 
mation, gy = ¢-, on the electrodes, giving the so-called ‘‘primary current 
distribution.’”’ Note, however, that the primary current distribution 
leads to unbounded current density at electrode edges or reentrant 
corners, while the secondary current distribution does not. The ex- 
ponential form of the electrode current relation, however, suggests 
that the primary and secondary current distributions differ only in a 
very small neighborhood of such a boundary singularity. 

Once ¢ has been determined, the electrode current densities can be 
calculated. Furthermore, the approximation does not depend on the 
specific electrolyte, but only on the existence of thin boundary layers 
across which the concentration changes are small. It thus applies to 
dilute, binary and concentrated, multicomponent electrolytes alike, 
although, of course, the various constants, e.g., the conductivities, may 
be vastly different in the two cases. 

To obtain a better approximation, we must include the effect of 
concentration variation across the boundary layers. Before doing so, 
however, let us consider a method of practical solution for ¢, the 
method of boundary integrals widely used for the numerical solution 
of potential problems. Since ¢ is harmonic, it satisfies Green’s third 
identity, which in the plane has the form 


are(P) = f | SE @GP, @ - o@ FP, @ | as(0, 


where [ is the curve bounding the electrolyte and G(P, Q) 
~ In (1/|PQ]), as the point P — Q, and is harmonic otherwise. When 
P lies on T, the identity becomes an integral equation for ¢, nonlinear | 
because d¢/dn is a nonlinear function of yg. Very efficient numerical 
methods for solving this equation have been developed. 


lil. THE BINARY ELECTROLYTE 


For the dilute, binary electrolyte, J. L. Blue? has extended the above 
boundary integral method to include the effect of concentration changes 
across the boundary layers. He derived a relation giving the concen- 
tration changes across the boundary layers in terms of integrals of 
the product of electrode current density and a kernel function whose 
form depends upon the local velocity distribution. These concentra- 
tion changes affect the electrode current density in two ways. First, 
as we have already noted, the ion exchange current density is a func- 
tion of concentration. Second, the potentials at the electrodes differ 
from those at the boundary layer edge, where the previous ¢ is actually 
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evaluated, by the so-called “concentration overpotentials’” Ag,, 
e = a, c. These overpotentials are due to the net current carried by 
ions diffusing at different rates through the boundary layers. To the 
present degree of approximation, they are linear functions of the 
concentration changes across the boundary layers. 

In a binary electrolyte, in which the two ion concentrations must be 
proportional, since the electrolyte is electrically neutral, Ag, is pro- 
portional to the integral giving the concentration change. This con- 
stant of proportionality, the constant giving the concentration to be 
inserted into the ion exchange current density, and the bulk conduc- 
tivity are the only constants containing material properties in the final 
relations. They thus characterize the electrolyte completely in this 
case and may thus be used in place of the individual diffusivities and 
mobilities. 

In the above, we did not qualify the term ‘‘binary electrolyte’ with 
the term ‘‘dilute.”? The discussion of the previous paragraph applies 
as well to concentrated binary electrolytes, provided that the boundary 
layer concentration changes are small. In this case, in the boundary 
layers, the equations are linearized around a state of constant bulk 
concentrations, so that the diffusivities, in general forming a 2 X 2 
matrix, and the mobilities, both functions of the concentrations, are 
constants, evaluated at the constant values of the bulk concentrations. 
Using the charge neutrality condition, the potential and one concen- 
tration are eliminated from the two boundary layer equations govern- 
ing the ion fluxes. The resulting equation has exactly the same form 
as that for a dilute, binary electrolyte, so that the solution has the 
same form, except that the proportionality constants have different 
values. 

In the following, we extend this result to concentrated, multicom- 
ponent electrolytes, the case of the ternary (n = 8) electrolyte being 
of particular interest. As above, we linearize the equations around a 
state of constant bulk concentrations, using the boundary layer 
approximation (not the assumption of excess supporting electrolyte). 
We eliminate the potential and one concentration to obtain a set of 
(n — 1) boundary layer equations for the concentrations, coupled by 
a reduced (n — 1) X (n — 1) diffusivity matrix. This system is solved 
by diagonalization, using the principal values and vectors of the re- 
duced matrix. Finally, these solutions are combined to calculate the 
electrode concentrations and concentration overpotentials. We find 
that both of these quantities are proportional to a single integral as 
before, so that the electrolyte can be characterized by specifying the 
two constants of proportionality, instead of n? diffusivities and n 
mobilities. 
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IV. THE COMPLETE BOUNDARY VALUE PROBLEM 


The steady flow of ions in a nonreacting, concentrated, n-component 
electrolyte is governed by the n conservation laws* 


V-N; = 0, (1) 


for 7 = 1, ---, n, where the ion flux N; is given in terms of the ion 
concentrations Ci, ---, C, and the potential ® by 


N; = CV = 3 D;;VC; = F2,U;C;V®. (2) 
j=1 


In these n relations, V is the circulation velocity of the electrolyte due 
to stirring, F is the Faraday constant, 2; the charge number of the 
ith ion, and the diffusivities D;; and mobilities U; are functions of 
Ci, «++, Cn.t We shall not impose any restrictions on the form of the 
D;;’s, although it seems reasonable to require, for example, that for 
1 # j, Di; —0, as C;—0. This makes the 7th ion flux N; vanish, 
when the ith ion is absent, and gives the correct form for a dilute 
electrolyte. With 

V-V =0, (3) 


as is appropriate for a liquid electrolyte, eqs. (1) and (2) yield the set 
of n nonlinear, second-order, partial differential equations 


v-( > D;;VC3 + FaUC.08) = V-VC;. (4) 

j=l 

The set of equations is completed by the charge neutrality condition 
x zC; = 0, (5) 
i=1 


which implies that the current density 


J= 00 FeN,; = —Kv@ — YF e,D;VC;, (6) 
i=1 i=] j=1 
with the conductivity 
K = > P2UC;, (7) 
i=l 
satisfies 
whe. (8) 


* Capital letters denote dimensional quantities; lower case letters denote dimen- 
sionless quantities. 

TIn Appendix C, we show how D,; and U; can be expressed in terms of the elec- 
trochemical potential. 
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We assume that there is no flux through the surface bounding the 
electrolyte, except that flux due to the production of cation 1 (e.g., 
Cutt) at the anode and deposition of the same cation at the cathode. 
Thus, for outward normal n, 


V-n=N;-n=0 (9) 
for? = 2, ---, n, while 
Ni-n = J-n/F2z (10) 
on the boundary, where 
—Ja, on the anode, 
Jun= J on the cathode, (11) 
0, otherwise. 


We shall assume that the electrode current densities determined by 
the surface reaction rates have the form 


Je = J(C1/C%) *Lexp(a-F|& — ®,|/RT) 
cas exp(—8.F|® = @,|/RT)], (12) 


for e = a, c, where the ion exchange current density J?(C1/C%)” is 
assumed to be a function only of the concentration of the cation 
participating in the surface reactions, R is the gas constant, 7’ the 
absolute temperature, ® the electrolyte potential, @, the electrode 
potential, and | — @,| the surface overpotential. These conditions 
complete the boundary value problem. 

Now let L, V, C, RT/F, D, D/RT, and FDC/L be some typical 
length, velocity, concentration, potential, diffusivity, mobility, and 
current density and set 


(x, y, 2) = (X, Y, Z)/L, 


C= C./C, 

g = F&/RT, 
v= V/V, 

j = LJ/FDC, 
di; = Dj;/D, 
UU = RTU;/D. 


The equations and boundary conditions then take the dimensionless 
form 


v-( 2 dijVe; + zusci¥ ¢) = Pev: Vc, (13) 
j=l 
for? = 1, ---, n, where the Peclet number Pe = VL/D, 
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> 2ic; = 0, (14) 
i=1 


v-(iV¢ os = >» eid) = 6, (15) 
where 
k= ~ ZU iC. (16) 
With v-n = 0 on the boundary, 
Ei dutelont eaedelom = PIG toe GL yy, 
where 
f= oe one = > aid; ;0c,/on 
_ a on the anode, 
=+— Je, on the cathode, (18) 
0, otherwise, 


and, for e = a, ¢, 
jeci*Lexp (ae| aa vel) = exp(—Be| gy — vel) J, (19) 
(LI(/FDC) (C/C). (20) 


Je 
ok 
Je 


V. THE BOUNDARY LAYER APPROXIMATION 


To obtain explicit results, we consider a plane, rectangular region 
with a cathode at Y = y = 0 and an anode at Y=L, y=1. We 
assume that the electrolyte is well stirred, with Pe > 1, so that the 
concentrations are nearly constant, except in thin boundary layers, 
of thickness of order eZ, « = (Pe)~}, balancing diffusion and convection 
near the electrodes. Far from the electrodes, we introduce the outer 
expansions 


é:(zx, Yy, €) a cf? oF O(e), 
a(x, Y; €) = e (a, y) - ep) (x, y) Sls 7) 


since eq. (13) implies that v- Vé; = O(e), where the c{’s are the con- 
stant bulk concentrations and ¢, g are harmonic functions; 1.e., 


92 + by = 2 + Sy = 0. (22) 
In the boundary layers, we set 


a= Cue, 7, €) _ cf? = eC}? (a, 9) 5 a) 
: G(x, Y, €) = 6 (x, 9) oF eG (x, 9) Se *y 


Ci 


Ye 


(21) 
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where the stretched variable 7 = y/e in the cathode boundary layer 
and 7 = (1 — y)/e in the anode boundary layer. The assumption that 
the electrode current densities are bounded has already been intro- 
duced by choosing the leading term in the expansion of ¢; to be con- 
stant. In the boundary layer eqs. (18) and (15) become 


(= ijt iz + eh) = (Wilie = Waliz) 
j= v 
= —eé ( > jC je + chides) , (28) 
j=l zx 


(Ee + » = eid: i) = —é¢ (te. + de 2% d.itie) ; (24) 
where we have introduced the reduced stream function 

U(x, 9, «-) = €*(z, ef), 
with the « and y components of v given by 


u= Py, v= —Wu, 


vanishing on the electrodes. The tildes on Pigs etc., mean di; 


= dij(G1, «++, En), etc. Note that the d’s, as, k, and J canst all be 
expanded around e = 0 to obtain the equations satisfied by ¢{?, 6, 
ge, +--+. The boundary conditions, for example at the cathode 


y = § = 0, have the form 


3 deta + eaten = [Ae OTD (25) 


0, for += 2,:+--,n, 
where 
je = jc(é1)”*Lexp(ae¢) — exp(—8-¢) ], (26) 


¢, and ¢ > 0 being evaluated at 7 = 0 and the cathode assumed to be 
grounded, so that y, = 0. Obviously 7, must also be expanded in 
powers of ¢. Finally, multiplying eq. (25) by z; and summing yields 


n 


a+ ae » 2:50 55 = Ee. (27) 
The set of relations from which the boundary layer equations are ob- 
tained is completed by matching Z;(z, 7, €) with é:(2, y, «) and ¢(z, 9, €) 
with 4(2, y, ¢) for small y and large g. Specifically, we use intermediate 
matching, setting 
y= ety*, Gg = ey", 
for fixed y* and 0 < a < 1, and considering the limit e — 0. Matching 


E(x, et 1y*, €) = cf + (a, et 1y*) 
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with 
éi(n, ery", 6) = of + O(8) 
trivially yields the boundary conditions 
Z(m, 0) = 
fori = 1, ---, n. Similarly, matching ¢ and g gives 
G(x, ©) = H(z, 0). 


The differential equation (24) and the pounds condition (27) 
imply that $5; o = ay = 0, so that | 


G(x, 7) = (a, 0). 
To a first approximation, the outer potential 4, as yet undetermined, 


simply penetrates the boundary layer. 


Vi. FIRST-ORDER MATCHING: THE SECONDARY CURRENT DISTRIBUTION 


Boundary conditions on ¢ are found by matching 


= 6 4+ eg@ 
with 
B= 6 + 6g, 


The zero-order inner potential ¢ has already been determined. The 
first-order potential ¢® satisfies the differential equation 


[> a9? + +d 2cdlPe |. = 0, 
L ] y 


t=1j7=1 


and the boundary condition 


bo oo os doze? = 7 
a) +l OD adipey Jey 


t=1 j=1 


at g = 0, where d{, k® are evaluated at cf, ---, c and j® at cf”, 
(x, 0) = (a, 0), so that 


jO(x) = Fle ]*{expla.e (x, 0)] — exp[—6.6 (2, 0) J}. 
Two integrations yield 


2(a, 1) = 0% (0) + [950@) — EY eager, n|/k, (28) 
t=1j=1 


where $“)(x) is an arbitrary function of x. The last term can be 
identified as the concentration overpotential. Since, for large g, ¢{? — 0, 


B~ G(x, 0) + [BO (a) + ety 72 (x) /k J, 
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for 7 = e*y* and «— 0. On the other hand, expanding in Taylor 
series gives 


o~ 6 (a, 0) + e*y* By (2, 0) + eG (a, 0) +... , 


for y = e*y*. First-order matching then yields the zero-order boundary 
condition 
By (#, 0) = j(2)/k 
and 
G(x) = (a, 0), 


another penetration condition. A similar analysis produces the 
boundary condition 
By (x, 1) = j(2)/k 


at the anode y = 1. Since j and 7 depend only on 4, we now 
have a complete boundary value problem for 4. To a first approxi- 
mation then, the potential ¢ satisfies Laplace’s equation 


Puz + Pyy = 0, (29) 


in the electrolyte, and the nonlinear boundary condition 


2a; on the anode, 
kag/dn = 4+— Je, on the cathode, (30) 
0, otherwise, 


where, for e = a, ¢, 
je = jeLe ]*Lexp(ae| gy — vel) a exp(—£.| g— vel), (31) 


yielding the so-called ‘“‘secondary current distribution.” 

Although we have derived this secondary current distribution 
specifically for small boundary layer thickness and moderate electrode 
current density, it probably is valid over a much wider range. The 
basic assumption is simply that the concentration variations are small. 
Thus, it is probably still valid when the flow is turbulent and the 
present simple description of the flow field does not suffice or at loci 
of boundary singularities (sharp edges, etc.) where the gradients are 
not necessarily normal to the boundary, as is assumed tacitly when 
one introduces the stretched variable 7 in the boundary layer approxi- 
mation. In any event, by calculating the first-order terms ¢{?, 3, 
we can obtain an error estimate or an improved approximation. 


Vil. CONCENTRATION VARIATIONS 


Although we have already derived boundary conditions on the ¢}””’s, 
we have not yet displayed the equations which they satisfy. With 
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constant leading terms c{° and ¥ expanded in the form 
V=Poteat-:- =Pfa) + o@+--:, 

where F(e) = Vur(2, 0) /2, g(x) = Voy (2, 0) /6, very €q. (23) becomes, 

to first order, 


> dE + (ki /2) 9 = POC — PO, (32) 


forz = 1, ---, n, where the 7th partial conductivity 
k= Zui , 


while the boundary condition (25) at 7 = 0 yields 


id gine? (0) Oa gO /2z,, for «+= 1, 
dd tj Cig =F (ki /z i) Gj | 0, for t = 2, ren (33) 


| 
— 


We can use the charge neutrality condition 


Ase = a a (25/2n)e9? 


and the relation 


-) _ 2¢ qe 
bg = jO- > % and DICH 


p=1 j= 


found from eq. (28), to eliminate ¢ and g“). We obtain the reduced 
boundary layer equations 


= digiy = Vet? — $Qaz, (34) 
forz = 1, ---, m — 1, and the boundary conditions at 7 = 0 
xo  |(l — t)9M/ei, for ¢ = 1, 
= digi = eee for 1=2,---,n—1, (35) 


where 4 = k{/k® is the ith transference number and the reduced 
diffusivity matrix [d;,;] has elements 


diy = UP — (2;/en)d2 — = (2p/z)[dp? — (25/2n)dmn ], (86) 


fort = 1,---,n —1, 7 = 1, ---,n — 1. The boundary value problem 
is completed by the conditions at infinity 


ci (x, ©) = 0, (37) 
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and the “‘initial condition”’ 
ef (0, 7) = 0. (38) 


As usual, the boundary layer approximation changes the equations 
from elliptic to parabolic type, requiring a condition at x = 0 in the 
time-like variable x. Equation (88) is equivalent to the assumption 
that the concentrations have their free-stream, bulk values near the 
leading edge x = 7 = 0. We now have a complete set of coupled 
equations and boundary conditions for the first-order boundary layer 
concentration variations ¢!”. We note in passing that the corresponding 
set for the second-order variations ¢{ would be much more com- 
plicated, containing not only ~“), as in the case of the dilute binary 
electrolyte, but also first-order terms d{?, etc., due to the concentration 
dependent diffusivities and mobilities. 


Vill. THE VECTOR CD-EQUATION AND ITS SOLUTION 


Equation (34) and boundary condition (35) can be written in the 
vector form 


[ae dgg = oP eM]. — YP]; (39) 
in g > 0, and 
Ca’ Le Js = (PO T9” (40) 
on g = 0, where [d’] is the (n — 1) X (nm — 1) matrix 
dis ae din—1 
ao ee ae 
din—11 ote din—1) (n—1) 
and [é ] and [7' ] the column vectors 
et? (1 — ti) /er 
reMwy=]o: f, (7To] = : 
e, —tna/2n—1 


By way of comparison, for the binary electrolyte (n = 2), one has 


~(1) 7(0)~ 7 (90)~(D 
tz = 0; CO -—YOG, (41) 
for g > 0, and 
~(1) . 
= TYID, (42) 


for g = 0, whered,, = 1,¢ = @, T® = TO (1 — t1)/z1. The solution 
of this problem is given in Appendix A, together with that of the 
corresponding axially symmetric problem, important for jet plating 
and for the rotating disk electrode. 
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We now introduce the eigenvalues d‘?) and the eigenvectors ['u‘”) ] 
of [d’ ], such that* 
Cd’ ][u] = d® [yu], 


for p = 1, ---, m — 1. If we set 
[zw] = = cu), 
p=1 


n—1 


[TO] = ru], 
p=1 


eqs. (89) and (40) are satisfied if 


(p) 7 (0) 710). ®) 
IM ee = bi oP — Wc; , (43) 
in g > 0, and 
-(0 
dre” ee 7(P) 7 y (44) 


on 7 = 0, forp = 1, ---,n—1. 
Now, in Appendix A we show that, if c(z, y) satisfies 


Cyy = Wylz — Wxly, 


with y = y’f(x), and 
c(0, y) = c(z, ~) = 0, 


then 
ae i ee OR Ga) det (45) 
where the kernel 
K@2)=-@'(f@Pae) /r@. (48) 
Thus, 
c()(x, 0) = [dir i * 5 (a) K (a, x")dx’, (47) 


with a factor (d‘?)3, found by replacing f by f/d‘ in K, and a factor 
(d‘?))— from the boundary conditions. Finally, to first order in «, 
the concentrations at the cathode are given by 


Z:(a,0) = of + ens f° J") Kea, @’)de’, (48) 
0 
for7 = 1, ---, m — 1, where 


n—1 
as = DY [dO 2 Ou, (49) 


p=l 


* We assume that such w’s and d’s exist. This must be verified in each particular 
case. 
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and we now write K, for K to indicate that the kernel depends on the 
velocity distribution near the cathode. A similar expression 


Gilt, 0) = ef + ens [5 (@")Kala, 2')de’, (50) 
0 
with x measured from the anode leading edge andz = 1, ---,n — 1, is 


valid for the concentrations at the anode. The concentration of the 
nth ion at cathode and anode is found from the charge neutrality 
relation 


in(m, 0) = — S (2i/en)8s(a 0). (51) 


With the Z{”’s determined, we are almost ready to calculate the first- 
order corrections to the electrode current densities. For example, at 
the cathode we have 

je = je(@1)“Lexp(ae¢) — exp(—6¢)], 
where 
2. = cf + ef, 
all quantities being evaluated at g = 0. With 
je I FIP Ho, 
we find 
j® = gle ]Lexp(a.6) — exp(—6.6) ], 
as before, and 
9? = FLAP P*(Lrae + a6] exp (aes) 
= [Ye cf» = Bop | exp(—6.¢)}, 
where 


0 = g(x, 0) — > 3 21d DE (x, 0) /k. (52) 


t=1 j=1 


We still must determine the first-order outer potential ¢“. It is the 
solution of the linear potential eee 


22 + Gy = 0, (53) 
in the electrolyte, 
qe); on. the anode, 
k9e®Y/an = 4 — 7, on the cathode, (54) 
0, otherwise, 


where the 7{?’s are linear functions of , with coefficients which 
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depend on cf, g, and é{?. This boundary condition is found from 
second-order matching of ¢ and ¢, just as the boundary condition 
on ¢®) was found from first-order matching. As is usual in boundary 
layer problems, the complexity of matching increases rapidly with 
order so this rather messy matching is relegated to Appendix B. Once 
@ has been determined, the electrode current densities can be calcu- 
lated to first order in e. 


IX. THE CONCENTRATION OVERPOTENTIALS 


The concentration variations ¢{? appear in the 7/”’s in two ways. 
First of all, there is the power law dependence of the ion exchange 
densities on the concentration ¢; of the ‘‘active’”’ cation (e.g., Cu**). 
This concentration is given by 


Pic eOdinae i * 7K da! (55) 
0 


in terms of the bulk concentration cf and the constant p1. 

The second concentration dependence is contained in the so-called 
“concentration overpotential,” the difference between the potential 
(x, 0) at the electrode and the potential g(x, 0) at the edge of the 
boundary layer. To first order in ¢, this concentration overpotential 
is given by 


Age = €(G@ — 6%) = 66 i "4K da’, (56) 
0 
where 


n n—l 
b= — DY aldip — (i/2n)dP Jui/k. (57) 
j=1 


i=l j= 


X. FORMAL PERTURBATION SOLUTION 


Although we shall eventually abandon this approach in favor of an 
equivalent scheme more convenient for numerical analysis, let us 
briefly summarize the successive steps by which the electrode current 
densities are calculated to first order in «. There are three such steps: 


(t) The nonlinear potential problem. First, we calculate the 
harmonic function ¢, satisfying nonlinear, mixed-boundary 
conditions. As is usual in boundary layer theory, we begin with 
a nonlinear problem, the only nonlinear problem in the se- 
quence. It yields first approximations to the electrode current 
densities, which may be adequate for many purposes. It also 
furnishes a boundary condition for the boundary layer con- 
centration variations. 
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(2) Concentration variations. The concentration variations é{”, or 
rather ¢{? and Ag., which are the only quantities involved in 
j@>, are calculated from eqs. (55) and (56), giving them in 
terms of a single integral of the product of electrode current 
density and given kernels. 

(zat) The linear potential problem. With ¢{, Ay. known, the har- 
monic function $"), satisfying linear, mixed-boundary condi- 
tions in ®, é{9, Ag, is calculated and, from it, a second ap- 
proximation to the electrode current densities, good to order e. 


XI. SOLUTION BY BOUNDARY INTEGRALS 


Instead of following the above rather cumbersome procedure, we 
can adopt the approach suggested by J. L. Blue,? forming a nonlinear, 
coupled set of integral equations whose solution, e.g., by successive 
approximation, yields the electrode current densities directly to order 
e. From this point of view, the formal perturbation scheme simply 
validates this method. For example, perturbation tells us that the 
integral solution of the boundary layer equation is good to order e 
and no more. 

We begin with Green’s third identity for a plane region: 


2re(P) = f Cia(Q)/KIGW, Qds(Q) — f Lie(Q/k] 
G(P, Qas(Q) — [ o(@ SP, QAO, 


where G has the fundamental singularity at P = Q, I. is the anode, 
I’, the cathode, and I the whole boundary, and we now denote 4 
+ee by g and the bulk conductivity k® by k. For each set of 
values of ja, je) ¢, and points P on the boundary, the solution of this 
integral equation gives ¢. Similarly, from the integrals 


Ci 


c+ ens I j Kd’, 
0 
Kova 2 / j-Kde’, 
; 0 


with e = a, c, we can calculate c; and Ag, from the given values of 
ja and j,. With ¢, ci, and Ay, now known on the electrodes, we can | 
calculate new values of j, and j, from the relations 


je = Je(cr) 7*Lexp (axe| e+ Ag — el) “— exp(—8.| g+ Ag — vel) J; 


for e = a, c, the anode potential ¢, > 0, and the cathode potential 
gy. = O (for a grounded cathode). Using these new values of electrode 
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current densities, compared with the old values, a successive approxi- 
mation scheme can be developed which converges rapidly. 

Note that this scheme is formally identical with that for the dilute, 
binary electrolyte. Furthermore, within the framework of this calcu- 
lation, the concentrated, n-component electrolyte can be characterized 
completely by giving its bulk conductivity k and the constants yw; and 
5. The values of the individual diffusivities d;; and mobilities u;, from 
which p; and 6 can be calculated, are of no importance. Only k, 1, 
and 6 enter the final set of equations. 


APPENDIX A 
Integral Solution of the Boundary Layer Equation 


' Let c(a, y) satisfy the boundary layer equation 


Cyy = WyCzr — Waly (58) 
inz>0,y > 0, with y = y’f(z) and 
c(0, y) = c(z, ©) = 0. (59) 


For use in the method of boundary integrals, we wish to express the 
boundary values of c in terms of an integral over the boundary values 
of cy, 1.e., we seek the kernel function K in the relation 


c(z, 0) = I ” oy (2", 0)K (a, 2")dz’. (60) 
0 
Following Levich? and Lighthill,‘ we make the change of variable 
g= [C¢@ Iz, 1 = = aI@T, 


the form of 7 being motivated by the observation that the right-hand 
side of eq. (58) vanishes when c is a function of y alone. We find that 
v(&, n) = c(2, y) satisfies the equation 


Vey = 2yv:, (61) 


if > 0, 1 > 0, and 
v(0, n) = v(f, ©) = 0. (62) 


Denote the boundary values of v, by g, so that 
v,(&, 0) = g(é). 


Then the Laplace transform 3, given by 
H(n, 8) = fo eet o(8, nds, 


BOUNDARY LAYER 819 


satisfies the ordinary differential equation 
vo’ — 2snd = 0, 

in0 < » < o, and the boundary conditions 

v’(0, s) = g(s),  (%,s) = 0, 
where g is the transform of g. Now suppose u(é, 7) is any solution of 
eq. (61), satisfying eq. (62). Then we have the identity between 
transforms 

O(n, 8s) = g(s)a(n, 8)/a'(0, s), 

giving 

5(0, s) = 0'(0, s)[u(O, s)/a’(O, s) ]. (63) 


A convenient choice for wu is the similar solution of eqs. (61) and (62), 
satisfying the condition 


u(é,0) = 1, 
so that a(0, s) = 1/s. In this case, u is a function of ¢ = 7/€ alone, 
given by 
= Ee exp (— 223/9) dz /(5) 7G) 
¢ 6 3)’ 
so that 


uy(é, 0) = —§-4/(8) TG) 
and, with I'(k)/s* the transform of #-4, 
aw (0, s) = —T(3)/(s)'P' (4)s?. 


Thus, 
a(0, s)/w’(0, s) = —(%)*P(3)/T(3)s*. 


Inversion of eq. (63), using the convolution theorem, then yields 
g 
v(&, 0) = f° og(2, Oe(e — Ede, 


where x, the inverse of @(0, s)/a’(0, s), is given by 


K(é) = — (¢)#e-4/1(3). 


Using the relations 


v(€, O)de’ = cy(2', O)dr’, 


[[Li@ Pa, 


Se 


we finally obtain eq. (60), with 


K(a, 2’) = — Gxes tia) * /r(3): (64) 
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K is of convolution type, i.e., a function of x — 2’ alone, only when 
f is constant, e.g., for Poiseuille flow near the wall of a channel or 
tube. For the laminar boundary layer flow on a semi-infinite flat plate, 


fas, 


so that 
K ~ (at — 2/84, 


A slight modification of the above yields the integral relation 
60) I " es(r’, O)L(r, r’)r’dr’, (65) 
0 


for the axially symmetric problem, where 

Coe = 1 (Pilr — Prez), 
inr > 0, 2z> 0, with y = 2?f(r) and 

c(0, z) = e(r, ©) = 0, 


where the kernel function 


Ler) = -(%) (fi tropa)” /r(Z)- (8) 


We merely replace the previously defined £ by 
e= |" C4) Foas. 


Equations (65) and (66) apply, for example, to the rotating disk elec- 
trode and to “jet plating,” where a jet of electrolyte is directed against 
an electrode. 


APPENDIX B 
Second-Order Potential Matching 


To obtain a boundary condition for the first-order outer potential 


@%), we must match 
B= GO + GO + FG 


with 
= GO + gO + 2g), 
where 
ge = g(x, 0), (67) 
EY = GO (a, 0) + GIP /KO — VY adiPeP/k. — (68) 


The second-order inner potential ¢®) satisfies the equation 


Je aP + Ema? + 3 XL ade + + dpe). = —k92, (69) 
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and the boundary condition at 7 = 0 
kB? + RWB? + > oy ad Pegg + apes) = J, (70) 
=a em 
where k, d:;, j- have been expanded in the forms 
=fO4 Z04.. ., 
a? + ed{P tones, 
jO + GM + --- 


Note that &®, d{P, 7{ are linear functions of the Z{’s; for example, 


di; 
Je 


k=hO ted (ak /ae,)2. 
p=1 


Two integrations give 


g® = p®) + GIP /kO ae Po ge /2 = PD > zids) EP) 


t=1 j=l 
+ f'[E & a€oapm — aye? — bo j9/0 | age, 
0 t=1j7=1 


where ¢@) is an arbitrary function of x and eq. (68) has been used to 
eliminate %. We wish to calculate ®@ in the limit 9 = e*—“y*, with 
0<a< 1, y* fixed, «— 0. The first three terms all contribute to 
the limit, while the fourth term vanishes, since ¢?) —-0, as 7 &. 
Finally, we must evaluate the limit of the integral. Consider the last 
term first. Since j© is a function of x only and 


RY) = = (0k /dc,)es”, 
we have n integrals like 
y 
[oe var. 

Changing variable gives 

ea—ly* y* 

[O° @ dag = et [PG ered > eT(0, 
0 0 


where I(e) +0, as «0. Similar limits are obtained for the other 
integrals with inteprands of the form ¢Z}7. Thus, 


B~ GO(H,0) + LH (a, 0) + em ty*jP/KO] + eG) 
+ et ly* 58> [k (0) _ (e%ly*)2 30) (x; 0)/2)], 


since 6 = ¢9 (zx, 0), to be matched with 


$22 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1976 


B~ B(x, 0) + ety* By" (a, 0) + (e*y*)* By (2, 0)/2 
+ [oO (x, 0) + ery* 3," (a, 0)] + B® (a, 0). 
Since ¢ is harmonic, g(x, 0) = — Sy(2, 0), and matching term by 
term gives the previous boundary condition on 9 
By (a, 0) = JP/b, 
the boundary condition on 4 
a (x, 0) = ae j/k, 
and the ‘‘penetration” condition 
6(x) = $2)(zx, 0). 
APPENDIX C 
The Electrochemical Potential 
In Section IV, we expressed the fluxes in terms of the gradients in 


the form 


N; = CV = 5 D;;VC; = F2z,C,U;V®, (71) 
j=l 


in effect defining the diffusivities D;; and mobilities U;, functions of 
the concentrations, by this relation. In electrochemistry, it is more 
usual to give the gradients of the electrochemical potentials, i.e., well- 
defined thermodynamic functions, in terms of the fluxes. Thus, New- 
man! gives the relations (for n + 1 species in our notation) 


ciVui = 2 Kii(V5 == V;,); (72) 
j= 
for 2 = 0, 1, ---, n, where the u,’s are the electrochemical potentials, 


the K;,’s are the interaction coefficients, functions of the concentra- 
tions, and V,’s are subject to certain constraints. Thus, 


3 CiVui = 0 
i=0 


and 


where V is the fluid velocity, p; = M;C; is the partial density for 
molecular weight M/;, and 
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These restrictions imply that 
ey (Ki; — K;i) = 0, 


for 7 = 0, 1, ---, n. In particular, Newman assumes the Stefan- 
Maxwell relations K;; = Kj, i.e., pairwise interaction between species. 
With these restrictions, eq. (72) has a unique solution for V; = V. 
Thus we find, following Newman, that 


V:-— Vo=- » LipCpVirp, 
= 
fori = 1, ---,n, where —[L,;] is the inverse of [M;; ], with 


= > Kix, for t= :F 
ae 
Ki;, for t¥ j. 


(Note that we may assume that K;; = 0, for 7 = 7.) Once V; — Vo is 
calculated, V; — V may be determined from the identities 


Ve ve = (ox/p) (Vi — Vo); 


VV; -— V=V;—- Vot Vo — V. 
With 


n 


Vip = DY (0np/0Cj) VC; + (du p/d) VE, 


j=0 


we finally find 


Do; = — »y 2X (px/p)LpCoC p(Oup/OC), 
= 


FzU) = — mu > (px/p)LipC p(dup/d®), 
=e 


Di; = > [Lis = ~ (ox/)La» | CC p(Ou,/9C)), 
Pats = ¥ [Lis - > (o1/)L| C'y(ptp/O8). 
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CHEMISTRY 


Amorphous Ferromagnetic Ag-X (X = Ni, Co, Gd) Alloys. J. J. Hauser, Phys. Rev. 
B, 12 (December 1975), pp. 5160-5165. The critical concentration for the appear- 
ance of ferromagnetism (41-at.% Ni) and the increase of 7’. with Ni concentration 
(~~9°K/at.%) for amorphous Ag-X (X = Ni, Co, Gd) films are quite close to the 
values reported for crystalline Ni-Cu alloys, suggesting again the validity of the 
virtual-bound-state model over the rigid-band model. 


Field Evaporation Experiments With a Magnetic Sector Atom-Probe FIM. T. Sakurai 
and KE. W. Miiller*, Surface Sci., 50 (May 1975), pp. 38-44. Field evaporation 
of Rh in the presence of *He and ‘He gases revealed that the formation of the helium 
compound (RbHe)?* is quite sensitive to He gas pressure; no helium compounds 
were observed below 5 X 10-7 Torr and all field-evaporated as helium compounds 
above 5 X 10-5 Torr at 78 K. * Pennsylvania State University. 


The Growth of LiNbO; Thin Films by Liquid Phase Epitaxial Techniques. A. A. 
Ballman, H. Brown, and P. K. Tien, J. Cryst. Grow., 29 (May 1975), pp. 289-295. 
A liquid phase epitaxial dipping process has been developed for the growth of thin 
films of LiNbO; on isostructural LiTaO; substrates. The process, employing a LizWO, 
flux, produces single crystalline layers from 1 to 10 um in thickness. Low-loss optical 
waveguides have been fabricated from these films which show a sharp step-change 
of index of refraction with the Li-TaQ3 substrate. 


MOS (Si-Gate) Compatibility of RF Diode and Triode Sputtering Processes. A. K- 
Sinha, J. Electrochem. Soc., 123 (January 1976), pp. 65-71. The radiation damage 
induced in the gate regions of CMOS capacitors by various thin-film sputtering 
processes has been evaluated. It is shown that trapped hydrogen can be used for 
annealing interface states in sealed MOS structures, which are impervious to hydro- 
gen during post-damage annealing. 


Photoreactions of Ketones With Amines. CIDNP Criteria for the Intermediacy of 
Aminoalkyl Radicals and Aminium Radical Ions. H. D. Roth and M. L. Manion, 
J. Amer. Chem. Soc., 97 (November 1975), pp. 6886-6888. The photoreactions 
of aromatic ketones with tertiary amines result in nuclear spin polarization and nmr 
line broadening. The intensity pattern of the polarization and the degree of line 
broadening depend upon the hyperfine coupling patterns of the paramagnetic inter- 
mediates and may distinguish between neutral aminoalkyl radicals and aminium 
radical ions resulting from net hydrogen abstraction and electron transfer, 
respectively. 


Preparation of BaTi;On by Solid-State Reaction. H. M. O’Bryan, Jr. and J. Thomson; 
Jr., J. Amer. Ceram. Soc., 58 (September—October 1975), p. 454. BaTi;Ou. forms 
as a reaction intermediate when a mixture of BaCO; + 4.5 TiO: is heated at 1000 
to 1100°C. Powder diffraction data are presented for BaTi;Ou. The conditions for 
forming this compound show similarities with those required to form the unstable 
BaTi,O; compound. 


Thermal Stability of PtSi Films on Polysilicon Layers. A. K. Sinha, 8. E. Haszko, 
and T. T. Sheng, J. Electrochem. Soc., 122 (December 1975), pp. 1714-1718. Films 
containing PtSi (2000 A) on poly-Si (1000 A) have been considered for application 
as first-level metallization in bipolar yst devices. The microstructure, resistivity, and 
stresses of these films were evaluated for various thermal treatments. No agglomera- 
tion was observed at up to 900°C; at 1000°C, a liquid-eutectic phase formed. 
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ELECTRICAL AND ELECTRONIC ENGINEERING 


Aberrations and Tolerances in a Double-Defiection Electron Beam Scanning System. 
M. G. R. Thomson, J. Vac. Sci. Technol., 12 (November—December 1975), pp. 1156- 
1159. A magnetic double-deflection system intended for use in an electron-beam lithog- 
raphy instrument has been theoretically analyzed. The first-order chromatic and 
third-order geometric aberrations are calculated so that the maximum field size can 
be found. The misalignment tolerance is estimated for different levels of required 
performance. 


An Annotated Bibliography of Microwave Circulators and Isolators 1968-1975. 
R. H. Knerr, IEEE Trans. Microw. Theory Tech., M7T'7-23 (October 1975), pp. 
818-825. A bibliography of microwave circulators and isolators from 1968 to 
1975 is presented (151 references). Some topics, which the author considered the 
most significant, include the eigenvalue approach to 3-port Y-junction circulator 
design, the peripheral mode isolator, electronic isolators, and wideband UHF circu- 
lators. The bibliography is subdivided into waveguide circulators, stripline circu- 
lators, lumped element circulators, other types of circulators, general papers appli- 
cable to all circulators, isolators, and books. 


Deep Level Spectroscopy, Low Temperature Defect Motion and Nonradiative 
Recombination in GaAs and GaP. C. H. Henry, J. Electron. Mater., 4 (October 1975), 
pp. 10387-1052). Recent development in junction capacitance measurements 
allow deep levels in semiconductors to be conveniently studied for the first time. 
The levels are often strongly coupled to the lattice which causes rapid nonradiative 
recombination. The released electronic energy causes violent vibrations near the 
defect which can promote low-temperature defect motion. 


Hopping Conductivity in Amorphous Carbon Films. J. J. Hauser, Solid State 
Commun., 17 (December 1975), pp. 1577-1580. The electrical resistivity of 
amorphous carbon films getter-sputtered at 95°K is well fitted between 300 and 20°K 
by the relation p = poexp[(T0/T7')?] with To ~7 X 10’ K. This behavior suggests 
a hopping conductivity very similar to that found in other amorphous semiconductors. 


Mass-Spectrometric Sampling of the Ionic and Neutral Species Present in Different 
Regions of an RF Discharge in Methane. M. J. Vasile and G. Smolinsky, Int. J. 
Mass. Spectrom. Jon Phys., 18 (October 1975), pp. 179-192. Ions originating in 
the sheath between the rr electrode and the plasma are characteristic of higher 
energy ionization processes than those observed from the plasma itself. Ion-molecule 
reactions observed in the rF sheath indicate a substantial proportion of encounters 
which occur at higher than thermal kinetic energy. 


A Per-Channel A/D Converter Having 15-Segment y-255 Companding. J. C. Candy, 
W. H. Ninke, and B. A. Wooley, IEKE Trans. Commun., COM-24 (January 1976), 
pp. 33-42. This encoder provides 8-bit compressed pcm and 13-bit uniform pcm 
at 8 kilowords per second. Also available is a 1-bit code that is useful for loop and 
switching applications. When suitably decoded the processed signal exhibits idle 
channel noise below 15 dBrnCo, harmonic distortion less than —40 dB, and gain 
tracking within +0.5 dB for circuit tolerances of 2%. The encoder is being integrated 
on two chips, and a decoder on a single chip. 


PHYSICS 


Adiabatic Following Model for Two-Photon Transitions: Nonlinear Mixing and 
Pulse Propagation. D. Grischkowsky,* M. M. T. Loy,” and P. F. Liao, Phys. Rev. A, 
12 (December 1975), pp. 2514-2533. The adiabatic following approximation is 
used to obtain solutions to a vector model describing a two-photon resonance. Explicit 
expressions are obtained for power-dependent nonlinear susceptibilities which are 
responsible for two-photon, resonantly enhanced, parametric processes. The results 
are also used to study the propagation of pulses nearly resonant with a two-photon 
transition. “IBM Thomas J. Watson Research Center. 
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Impact Ionization Rates for Electrons and Holes in GaAs;_,Sb, Alloys. T. P. Pearsall, 
R. E. Nahory, and M. A. Pollack, Appl. Phys. Lett., 28 (April 1, 1976), pp. 403-405. 
We have measured the impact ionization rates for electrons and holes in GaAs;_,Sbz 
for alloys with z = 0.05, 0.10, and 0.12. The ionization rates for electrons (a) and 
holes (8) are different and a function of Sb content. For GaAso.ssSbo.12, a/8 = 2.5, 
while for GaAs, a/8 = 0.25. 


Ion-Ion Correlations and Diffusion in Beta Alumina. D. B. McWhan et al., Phys. 
Rev. Lett., 35 (October 6, 1975), pp. 953-956 and 36 (February 9, 1976), p. 344(e). 
A comparison of the correlation length for short-range order (5, 6, 11, 18 A), the 
attempt frequency for diffusion (56, 28, 80, 85 cm~), and the previously reported 
values of tracer diffusion coefficient for Na, Ag, K, and Rb, respectively, suggests 
that ion-ion correlations inhibit diffusion in the 6-aluminas. 


Millimeter and Submillimeter Wave Laser Action in Symmetric Top Molecules 
Optically Pumped Via Perpendicular Absorption Bands. T. Y. Chang and J. D. 
McGee, IEEE J. Quantum Electron., Q#-12 (Jaunary 1976), pp. 62-65. Ninety- 
nine new far infrared laser lines from 227 um to 1.965 mm have been observed in 
CH;CN, CH;CCH, CH;Cl, CH3Br, and CHsI by optically pumping these gases with 
CO, laser pulses of 150-us duration. 


Motion of pn Junctions in CuInSe>. B. Tell, Sigurd Wagner, and P. M. Bridenbaugh, 
Appl. Phys. Lett., 28 (April 15, 1976), pp. 454-455. We report the first pn junc- 
tion delineation and diffusion study in a ternary chalcopyrite-type semiconductor. 
Pn junctions were formed in Zn or Cd plated p-CulnSez by five-minute anneals at 
200 to 450°C. By angle lapping and staining techniques, junction depths 2; varying 
from ~1 pm to ~130 um were determined, from which the interdiffusion coefficient 
derived is D(cm?/s) = 164 exp[—1.19(eV)/kT]. The large preexponential term in- 
dicates concentrations or mobilities of point defects substantially above those of 
the related II-VI compounds. 


Novel Technique for Measuring Nitrogen Profiles in GaP: N. J. Shah, R. F. Leheny, 
and P. D. Dapkus, J. Appl. Phys., 46 (December 1975), pp. 5244-5247. A high- 
resolution technique for determining the nitrogen-doping profiles in GaP:N using a 
laser-pumped dye cell as a broad-band fluorescence source is described. Results ob- 
tained suggest that, under some LPE growth conditions, vapor-liquid equilibration 
may be slow and, as a result, affect the nitrogen incorporated into the solid. The 
implications of these results for high-efficiency LED’s are also discussed. 


Observation of a New Nonlinear Photoelectric Effect Using Optical Levitation. 
A. Ashkin and J. M. Dziedzic, Phys. Rev. Lett., 36 (February 2, 1976), pp. 267-270. 
We observe a new three-photon nonlinear photoelectric effect in glass with visible 
cw laser light using optical levitation of transparent particles. Electrons are pumped 
into the conduction band of glass by the absorption of two photons and subsequently 
ejected in a single photon step. 


Polarization Rotation Induced by Resonant Two-Photon Dispersion. P. F. Liao and 
G. C. Bjorklund, Phys. Rev. Lett., 36 (March 15, 1976), pp. 584-587. A polari- 
zation rotator is demonstrated which utilizes the dispersion associated with 35 to 58 
two-photon transitions in sodium vapor. A linearly polarized beam at » is rotated 
by a circularly polarized beam at v2, and »: + v2 is near the two-photon transition 
frequency. When combined with a polarizer, a rapid optical shutter is obtained. 


Preparation and Properties of InP/CdS and CulnSe:/CdS Solar Cells. J. L. Shay 
et al., Proc. 11th IEEE Photovoltaic Specialists Conf., 1975 (Scottsdale, Arizona, 
May 6-8), pp. 503-507. Described are the preparation and properties of the 
recently reported InP/CdS and CulnSe2/CdS single-crystal solar cells. Solar power 
conversion efficiencies of 12.5 and 12% have been observed. The ultimate efficiency 
achievable with InP/CdS using state-of-the-art liquid encapsulated Czochralski 
p-type InP is calculated to be 17.2% for AM2 conditions and 14.0% for AMO 
conditions. 
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Push-rod AC Susceptibility Apparatus. J. J. Hauser and C. M. Antosh, Rev. Sci. 
Instrum., 47 (January 1976), pp. 156-157. A push-rod ac susceptibility apparatus 
which permits the measurement of ultrathin films of superconductors and ferro- 
magnets as a function of temperature and magnetic field is described. Films prepared 
at 77° K can be mounted under liquid nitrogen and measured without warmup. 


Water and Its Relation to Broken Bond Defects in Fused Silica. R. H. Stolen and 
G. E. Walrafen,* J. Chem. Phys. 64 (March 15, 1976), pp. 2623-2631. Vibrational 
modes associated with OH in fused silica were studied by Raman scattering. In- 
tensities of lines at 604 and 490 cm decrease with increasing OH and increase with 
fictive temperature, suggesting that OH is trapped at defect sites. A weak 970 cm™ 
line was seen in fibers and assigned to the Si-(OH) stretch. * Dept. of Chem., 
Howard University. 
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